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CHAPTER  1 


INTRQPVCTIQJS: 

The  objective  of  this  project  is  to  develop  an  electromagnetic 
wave  propagation  technique  for  the  assessment  of  the  subsurface 
condition  of  reinforced  concrete  structures  using  ground  penetrating 
radars.  Concrete  bridge  decks  have  been  studied  in  particular 
because  they  have  been  the  focus  of  other  ongoing  research  at  MIT 
w'hich  have  provided  actual  field  data.  In  addition,  concrete  bridge 
decks  incorporate  many  characteristic  features  of  other  reinforced 
concrete  structures,  and  thus  the  specific  results  obtained  from  this 
study  can  be  generalized  without  much  difficulty  to  other  related 
applications. 

Ground  penetrating  radar  is  the  electromagnetic  analog  of 
sonar  and  ultrasonic  pulse-echo  methods.  The  radar  antenna  can  be 
either  in  contact  with  the  bridge  deck  slab  or  10  to  30  centimeters 
above  it.  A  short  (0.5  to  1.5  ns.)  pulse  of  electromagnetic  energy  is 
transmitted  into  the  slab,  which  then  gets  reflected  from  objects  or 
interfaces  which  represent  discontinuities  in  electrical  properties 
(Maser,  1989).  In  a  bridge  deck,  such  discontinuities  include  the 
asphalt/concrete  interface,  top  and  bottom  rebar  grid,  the  bottom  of 
the  slab,  and  abrupt  variations  in  concrete  properties.  The  reflected 
echo  signal  is  then  picked  up  by  a  receiving  antenna,  which  is  usually 
same  as  -  or  located  adjacent  to  -  the  transmitting  antenna  (Maser, 
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1989).  Another  typical  characteristic  of  commercial  antennas  is  that 
they  are  linearly  polarized.  This  means  that  they  will  tend  to  show  a 
much  stronger  response  to  rebars,  cylinders,  or  other  subsurface 
anomalies  oriented  parallel  to  the  polarization  direction  of  the 
electric  field.  The  time  history  of  these  reflected  signals  is  referred 
to  as  the  "waveform". 

The  strength  of  the  peaks  in  this  waveform  (which  correspond 
to  the  subsurface  interfaces  and  anomalies)  depend  on  the  reflection, 
transmission  and  attenuation  characteristics  of  the  medium.  The 
reflection  and  transmission  coefficient  at  an  interface  depends  on  the 
velocity  of  the  layers  on  either  side,  which  in  turn  depends  on  their 
dielectric  properties.  The  relative  position  of  the  peaks  in  the 
waveform  also  depends  on  these  velocities.  Thus,  velocity  and 
attenuation  of  electromagnetic  waves  in  the  medium  are  the  two 
main  parameters  affecting  the  output  waveform,  and  these 
parameters  are  related  to  the  complex  dielectric  permittivities  of  the 
various  underlying  layers. 

During  the  first  year  of  this  project  (January  to  December, 
1987)  the  focus  of  research  was  on  modelling  the  reflection  of 
electromagnetic  waves  by  metallic  cylindrical  objects  buried  in 
concrete  and  on  the  synthesis  of  radar  waveforms  from  concrete 
bridge  decks  (Maser  and  Halabe,1988).  Experiments  were  carried 
out  to  observe  reflections  from  rebars  and  cylinders  suspended  in  air 
and  buried  in  sand.  These  reflections  were  then  modelled  using  a 
Scattering  Attenuation  Function  (SAF)  which  was  derived  from 
principles  in  geometric  optics.  This  model  for  rebars  was  then 
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incorporated  in  a  computer  pro;*.ram  for  waveform  synthesis.  The 
dielectric  properties  of  concrete  were  modelled  using  a  simple  three- 
phase  mixture  formula  and  radar  waveforms  were  synthesized  for 
normal  concrete,  low-moisture  deteriorated  concrete  and  for  high- 
moisture  deteriorated  concrete.  A  comparison  of  these  synthesized 
waveforms  with  radar  waveforms  obtained  from  actual  bridge  decks 
showed  that  the  field  waveforms  incurred  considerably  higher 
attenuation  than  that  predicted  by  the  synthesis  model.  Thus,  there 
was  a  need  for  improved  understanding  of  the  velocity  and 
attenuation  characteristics  of  electromagnetic  waves  in  reinforced 
concrete. 

The  objective  in  this  year’s  research  has  been  to  characterize 
the  velocity  and  attenuation  of  electromagnetic  waves  in  reinforced 
concrete  as  a  function  of  frequency,  temperature,  moisture  content, 
chloride  content  and  concrete  mix  constituents.  This  objective  has 
been  pursued  by  adaptation  of  theoretical  models  which  have  been 
developed  and  verified  for  partially  saturated  rocks.  In  these  models, 
the  electromagnetic  properties  of  the  medium  has  been  modelled  by 
aggregating  the  properties  of  its  most  significant  constituents.  The 
research  described  herein  has  adapted  these  aggregating  techniques 
to  reinforced  concrete,  whose  significant  constituents  are  coarse 
aggregate,  fine  aggregate,  cement  paste,  air,  water,  salt  and  steel 
rods.  Although  the  volumetric  content  of  water  in  a  concrete  mixture 
is  usually  small  compared  to  the  volume  of  aggregates,  it  has  the 
most  significant  contribution  towards  the  dielectric  permittivity  of 
the  overall  mixture  because  of  its  high  complex  dielectric 
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permittivity.  This  electromagnetic  property  of  water  is,  in  turn, 
highly  sensitive  to  the  amount  of  dissolved  salts.  Hence,  emphasis  in 
this  study  has  been  placed  on  modelling  of  dielectric  properties  of 
saline  water. 

The  following  chapters  will  describe  the  various  models  that 
have  been  considered,  and  the  numerical  results  obtained. 
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CHAPTER  2 

MODELLING  OF  DIELECTRIC  PROPERTIES 

2.1  Introduction 

A  dielectric  material  is  a  material  which  conducts 
electromagnetic  waves,  and  is  commonly  referred  to  as  an 
'insulator'.  The  transition  from  insulator  (which  is  characterized 
by  low  electric  conductivity)  to  conductor  (which  has  high 
conductivity)  is  a  gradual  one  without  any  clear  cut-off  point,  and 
many  existing  materials  exhibit  both  properties.  If  the 
conductivity  of  a  material  is  high,  it  attenuates  electromagnetic 
waves  to  a  larger  extent,  thus  resulting  in  a  lower  penetration 
depth.  The  two  significant  properties  of  a  dielectric  material  are 
the  real  and  imaginary  parts  of  its  complex  dielectric  permittivity. 
The  real  part  is  associated  wih  the  phase  velocity  and  the 
imaginary  part  signifies  the  conductivity  or  attenuation  of 
electromagnetic  waves  in  the  medium.  The  dielectric  permittivity 
is  often  referred  to  as  the  dielectric  'constant'.  This  term  is  a 
misleading  one,  since  this  property  of  a  dielectric  material  varies 
with  several  factors.  Bell,  et.  al.  (1963)  have  provided  some 
insight  into  these  variations  through  underlying  molecular  and 
atomic  mechanisms  associated  with  dielectric  constant  and 
conductivity  of  such  materials. 

A  dielectric  material  increases  the  storage  capacity  of  a 
capacitor  by  neutralizing  charges  at  the  electrode  surfaces.  This 
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neutralization  can  be  imagined  to  be  the  result  of  the  orientation 
or  creation  of  dipoles  opposing  the  applied  field.  Such  a 
polarization  is  proportional  to  the  polarizability  of  the  dielectric 
material;  that  is,  the  ease  with  which  it  can  be  polarized. 
Polarizability  is  defined  as  the  average  induced  polarization  per 
unit  field  strength.  The  greater  the  polarizability  of  a  material,  the 
greater  will  be  its  dielectric  constant. 

There  are  four  recognized  mechanisms  of  polarization;  (a) 
electronic,  (b)  atomic,  (c)  orientation,  and  (d)  space  charge  or 
interfacial  mechanism.  The  first  three  polarization  mechanisms 
are  forms  of  the  dipole  polarization  previously  mentioned.  They 
result  from  the  presence  of  permanent  dipoles  or  from  dipoles 
induced  by  external  fields.  The  fourth  type  •  space  charge 
polarization  -  results  from  charge  carriers  that  are  usually  present 
and  more  or  less  free  to  move  through  the  dielectric.  When  such 
carriers  are  impeded  in  their  movement  and  become  trapped  in 
the  material  or  at  interfaces,  space  charge  concentrations  result  in 
a  field  distortion  that  increases  the  capacitance. 

The  above  discussions  indicate  that  the  dielectric  constant  of 
a  material  is  not  such  a  constant,  but  is  a  function  of  the 
polarizability  which  is  in  turn  a  function  of  frequency, 
temperature,  local  fields,  applied  field  strength,  availability  and 
freedom  of  charge  carriers  within  the  dielectric,  and  local  field 
distortions. 

The  dielectric  conductivity  (defined  later  in  Eq.  7)  is  also 
affected  by  the  same  factors.  This  is  because  the  dielectric 
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conductivity  is  a  function  not  only  of  ohmic  conductivity  but  also 
of  the  power  consumed  in  polarizing  the  material. 

In  addition  to  the  'intrinsic  loss'  caused  by  the  conduction 
process,  there  could  also  be  'scattering  loss'  due  to  the  presence  of 
inhomogeneities  within  the  medium.  Since  the  aggregate  particle 
size  is  much  smaller  than  the  electromagnetic  wavelength 
considered  here  (X  =  10  cm  in  concrete),  the  inhomogeneities  in 
concrete  are  not  likely  to  cause  significant  scattering  losses. 
However,  this  is  not  true  for  inhomogeneities  with  large  dimension 
(like  delaminations)  or  with  large  conductivity  (like  rebars).  These 
inhomogeneities  will  be  modelled  separately  and  their  effect  will 
be  taken  into  account  for  the  purpose  of  waveform  synthesis. 

2.2  Models  for  Dielectric  Properties  of  Saline  Water 

Moisture  and  chloride  are  the  two  constituents  which  have  the 
most  prominent  effect  on  the  dielectric  properties  of  concrete. 
Although  the  volumetric  content  of  water  in  a  concrete  mixture  is 
small,  it  has  a  very  important  bearing  on  the  velocity  and 
attenuation  of  electromagnetic  waves  in  concrete  because  of  its  high 
complex  dielectric  permittivity.  Repeated  application  of  de-icing 
salts  over  concrete  bridge  decks  in  winter  months  along  with  the 
cracking  of  the  concrete  due  to  freeze-thaw  cycles  cause  the  salt  to 
infiltrate  into  the  bridge  decks.  This  salt  is  commonly  referred  to  as 
free  chlorides  because  it  can  readily  go  into  solution  (as  opposed  to 
chloride  components  which  are  chemically  bound  in  the  cement 


-  8  - 


paste).  Only  the  chlorides  which  are  actually  present  in  solution 
form  at  any  given  time  affects  the  dielectric  properties  of  water  and 
of  the  concrete  mixture  as  a  whole.  Presence  of  dissolved  chlorides 
reduces  the  real  part  of  the  complex  dielectric  permittivity  of  water 
(which  affects  the  radar  velocity),  and  greatly  increases  the 
imaginary  part  (which  affects  the  attenuation  of  radar  waves).  The 
latter  effect  is  due  to  the  increase  in  ohmic  conductivity  of  water.  In 
addition,  temperature  of  water  affects  its  conductivity,  and  is  thus 
another  factor  affecting  its  dielectric  properties.  A  model  for 
predicting  these  dielectric  properties  of  water  is  presented  in  the 
next  paragraph. 

Various  researchers  (e.g.,  Dorsey,  1940;  Saxton,  1952b;  Stogryn, 
1971)  have  formulated  empirical  models  to  fit  the  observed 
experimental  results  for  dielectric  properties  of  aqueous  sodium 
chloride  (NaCl)  solutions.  Some  of  the  predicted  values  have  been 
tabulated  by  Dorsey  (1940)  and  Saxton  (1952b).  Earlier  models 
presented  by  Dorsey  (1940)  prescribe  different  equations  for 
different  temperatures,  and  thus  are  difficult  to  use.  The  model  by 
Saxton  (1952b)  is  also  difficult  to  use  because  of  its  complicated 
form.  Recent  work  in  this  area  have  removed  these  problems  and 
have  improved  the  accuracy  of  prediction  (Ulaby  et  al,  1986).  Ulaby 
et  al  (1986)  have  presented  models  developed  recently  which  are 
discussed  below.  These  are  empirical  models  obtained  by  a  non¬ 
linear  regression  analysis  to  fit  the  experimental  data  compiled  by 
several  researchers. 
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A  model  for  pure  and  saline  water  wa;.  developed  by  Stogryn 
(1971)  and  later  modified  by  Klein  and  Swift  (1977).  Since  tap  water 
or  water  in  the  pores  of  concrete  is  never  pure,  the  pure  water  model 
is  not  of  much  interest  here.  The  salinity  (Sg^)  of  a  salt  solution  is 

defined  as  the  total  mass  of  solid  salt  in  grams  dissolved  in  one 
kilogram  of  the  aqueous  solution.  Thus,  Sjw  is  expressed  in  parts  per 
thousand  (ppt)  by  weight.  The  normality  of  saline  water  (Njvv)  is 
related  to  its  salinity  by  the  following  expression  (Klein,  1977): 

Nsw  =  A  Ssw  11-707  X  10*2  +  1.205  x  10*5  Sgw  +  4.058  x  lO’^  Sgw^] 

(1) 

where  A=1  for  sodium  chloride  (NaCl)  solution,  and  A=0.9141  for  sea 
water,  which  includes  other  salts  as  well.  The  above  equation  is 
valid  for  Sjw  <  260  ppt. 

The  frequency  dependence  of  the  complex  dielectric  constant 
(^sw)  of  saline  water  is  given  by  the  well-known  Debye  equation 

(Ulaby,  1986): 


£sw  “  Eswoo  + 


£swo  ~  £swgo 

1  -  i  27Cflsw 


+ 


i  Pi 
27C£of 


(2) 


where, 

^swo  =  static  dielectric  constant  of  saline  water  (dimensionless) 
£ 

swoo  =  very  high  frequency  (or  optical)  limit  of 
(dimensionless) 

X 

sw  =  relaxation  time  of  saline  water  (seconds) 


10  - 


%  • 


f  =  electromagnetic  frequency  (Hz) 

- 1  2 

Eo  =  dielectric  permittivity  of  free  space  =  8.854  x  10 
farad/m 

c\  =  ionic  (i.e.,  ohmic  or  dc)  conductivity  of  the  aqueous 
saline  solution  (mho/m) 


Upon  rationalizing  Eq.  (2),  can  be  expressed  in  terms  of  real  and 
imaginary  part?  as; 

^sw  “  ^  sw  ^  sw  (3) 

where  the  real  and  imaginary  parts  are  given  by: 


Esw  ”  Eswoo  + 


Eswq  ~  £swoo 

1  +  (27CfTsw)2 


(4) 


and 


__  27CfXsw(£swO  "  Eswoe)  J 
Csw  ^  ^ ^ 


Oi 


1  +  (27CfXsw)2  27C£of 


(5) 


The  imaginary  part  (e''sw)  is  also  termed  as  the  'loss  factor'.  For  pure 
water  oj  is  zero  and  hence  the  second  term  in  Eq.  (5)  vanishes.  The 
loss  tangent  (also  called  'dissipation  factor’)  in  w'ater  is  defined  as: 


uui  5w  —  “1*^  (6) 

Esw 

The  dielectric  conductivity  of  water  is  defined  as: 


=  (D  Es>v  Eo  (7) 

where  to  is  the  angular  frequency  i-lnf).  This  conductivity  is  a 

frequency  dependent  quantity  which  can  be  caused  by  the 
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movement  of  charged  bodies  such  as  ions,  the  orientation  of  dipoles, 
or  any  other  process  of  energy  dissipation  (Bell  et  al,  1963),  It  can 
be  seen  from  Eqs.  (5)  and  (7)  that  in  the  d.c,  limit  (cd  0)  the 
dielectric  conductivity  (Ow)  is  nothing  but  the  ionic  conductivity  (Oj) 
of  water. 

According  to  Stogryn  (1971)  there  is  no  evidence  to  indicate 
that  Cg  depends  on  salinity,  and  a  value  of  Eg  =  4,9  provides  a 

good  fit  to  the  available  data.  The  dependence  of  Egvvo  ^sw  (PPO 

and  temperature,  T  (^C)  can  be  obtained  using  the  following 

empirical  expressions  developed  by  Klein  and  Swift  (1977): 

^swo  ^sw)  ~  ^swo  ^sw)  (^) 


where. 


eswo(T,0)  =  87.134  -  1.949  x  IQ-^  T  -  1.276  x  lO'^T^ 

+  2.491  X  10*4  (9) 

a(T,  Sgw)  =  1.0  +  1.613  x  10*5  T  Sgw  -  3.656  x  10-3  Sg^v 

+  3.210  X  10-5  85^2  .  4.232  x  10-7  (10) 


The  above  equations  are  valid  for  salinities  in  the  range  4  <  Sgw  <  35 
ppt.  The  form  used  in  Eq.  (8)  may  also  be  used  to  define  the 
relaxation  time. 


'^sw  (T,  Ssw)  “  “^sw  (TtO)  .  b  (T,  Sgw) 


(11) 
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Here  Xgvv  is  the  relaxation  time  for  pure  water  and  can  be 

obtained  using  the  following  expression  presented  by  Stogryn 
(1971): 

2n  Tsw  (T,0)  =  1.1109  X  10-JO  -  3.824  x  lO'J^T  +  6.938  x  lO-l^T^ 

-  5.096  X  10* (J2) 

A  related  term  used  in  the  literature  is  the  relaxation  frequency  (f^), 
where. 


fr=  (271  X  )-  J 
f  sw 


(13) 


It  can  be  seen  from  Eq.  (5)  that  the  first  component  of  e"sv^'  has  i:s 
maximum  value  at  f  =  fj-.  The  expression  for  b(T,  Sjw)  '^as  first 
presented  by  Stogryn  (1971)  and  later  modified  by  Klein  and  Swift 
(1977)  through  the  use  of  Eq.  (1)  to  yield: 

b(T,  Ssvv)  =  1.0  +  2.282  x  10*5  T  Sgw  -  7.638  x  10-4  Sgw 

-  7.760  X  10'^  Ssw^+  1.105  x  10'^  Ssw^  (14) 

This  expression  is  based  on  data  for  0  ^  T  ^  40^C  and  0  ^  Ngw  ^  3  (or 
approximately,  0  ^  Sjw  ^  157  ppt  for  an  NaCl  solution). 
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The  only  parameter  that  now  remains  is  the  ionic  conductivity, 

o  .  For  sea  water,  o  was  derived  by  Weyl  (see  Stogryn,  1971)  and 
i  i 

later  modified  by  Stogryn  (1971)  to  the  following  form: 

o  (T,  S  )  =  a.  (25,  S  )  e'*?  (15. 

1  sw  1  sw 

where  g.  (25,  S^^)  is  the  ionic  (i.e.,  ohmic  or  dc)  conductivity  of 
saline  water  at  25  °C  and  is  given  by: 

Oj  (25,  10.18252  -  1.4619  x  lO'^  +  2.093  x  lO'^ 

-  1.282  X  lO’’^  S  (16) 

SW 

The  function  <p  in  Eq.  (15)  depends  on  Sjvv  A  =  25  -  T, 

(p  =  A  [2.033  X  10-2  +  1.266  x  10-4  A  +  2.464  x  10-6  a2 

-  S^^(1.849  X  10'^  -  2.551  x  10'^  A  +  2.551  x  10'^  A^)]  (17) 

s  w 

The  above  expressions  are  valid  in  the  range  0  ^  -  40  ppt. 

Different  expressions  presented  in  this  section  are  valid  over 
different  salinity  ranges.  The  combination  of  these  equations  is  valid 
in  the  range  4  ^  Sg^v  -  35  ppt  and  0  ^  T  S  40  ^C.  The  above  model 
can  also  be  used  for  pure  water  (Sgw  =  0)  with  one  small  variation. 
For  pure  water,  Eq.  (9)  for  Eg  wo  1°  6e  replaced  by  the 

following  expression  given  by  Klein  and  Swift  (1977). 

Eswo  (pure  water)  =  88.045  -  0.4147T  +  6.295  x  10-4t2 


+  1.075  X  lO'^T^ 


(18) 


Loss  Factor  of  Mater 
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Frequency  (6H2) 

(D) 

Figure  1  (a)  Real  and  (b)  Imaginary,  Parts  of  the 

Complex  Dielectric  Constant  of  Water  (from 
Ulaby  et  al ,  1986) 
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The  spectral  variations  and  dielectric  relaxation  phenomena  for  e’ 

s 

and  £"sw  are  shown  in  Figure  1  for  pure  water  and  saline  water  with 
Ssw  =  32.54  ppt,  which  is  an  average  value  for  sea  water  (Ulaby, 
1986).  The  relaxation  frequency  ((q)  of  pure  water  is  nearly  9  GHz 
at  O^C  and  17  GHz  at  20°  C.  It  can  be  seen  from  Figure  1  that 
presence  of  dissolved  chlorides  has  a  very  prominent  effect  on  the 
loss  factor  of  water,  particularly  in  the  'radar'  range  (0.1  to  10.0  GHz). 

The  above  model  (Sgvv-  <  35  ppt)  is  adequate  for  concrete 

exposed  to  a  marine  environment.  However,  a  significantly  higher 
salinity  (of  the  order  of  140  ppt)  has  been  reported  for  concrete 
bridge  decks  (see  compulations  in  the  chapter  on  'numerical 
studies').  The  chloride  content  at  saturation  is  357  ppt  for  cold  water 
and  391.2  ppt  for  hot  water  (Weasi,  1986).  Water  with  such  high 
salinity  is  commonly  referred  to  as  'brine'.  A  model  for  dielectric 
permittivity  of  brine  was  developed  by  Stogryn  (1971)  and  is 
presented  below.  This  model  has  been  developed  to  fit  the  observed 
experimental  data  in  the  range  0  ^  Ngw  -  3  (i.e.,  0  ^  Sgw  ^  157  ppt 
for  NaCl  solutions)  and  0  ^  T  ^  40°C.  However,  for  Sgv^r  ^  35  ppt,  the 
equations  presented  earlier  in  this  section  give  more  accurate 
predictions  (Klein,  1977). 

In  terms  of  Siogryn's  (1971)  formulation  for  brine: 

(19) 


^swo  ("!'>  l^sw)  -  ^swo  •  ^1  (l^sw) 
^sw  l''sw)  ~  “^sw  ('1'’®)  •  1^1  l^sw) 


(20) 


-  16  - 


("Ti  ^sw)  “  ^sw)  •  ^sw)  (21 ) 

The  above  functions  are  given  by  the  following  expressions: 

^swo  (T,0)  =  *^swo  (pure  water)  given  by  Eq.  (18), 
ai(Nsw)  =  1.0  -  0.255  Ngw  +  5.15  x  10-2  Ngvv^ 

-  6.89  X  10-3  Ksw^  (22) 

“^sw  (T,0)  is  defined  by  Eq,  (12). 

t>l  (T,  Nsw)  =  1.0  +  0.146  X  10-2  t  N^vv  -  4.896  x  10-2 

-  2.97  X  10"^  +  5.64  x  10*^  Nsw^ 

Combining  Eqs.  (1)  and  (16),  with  A=1  for  NaCl  solution,  we  obtain: 

Ci  (25,  Nsw)  =  Nsw  (10.39  -  2.378  Kjw  +  0.683  Ngw^  -  0.135  Ngw^ 

+  1,01  X  10'^  Nsw"^)  (24) 

The  last  parameter  cj  depends  on  Nsvv  and  A  =  25  -  T, 

Cj  (A,  Nj^)  =  1.0  -  1.96  X  10*2a  +  8.08  x  lO'^  a2 

-  A  [3.02  X  10-5  +  3.92  x  10'^  A 

+  Nsw  (1.72  X  10-5  -  6.58  x  10-6  A)]  (25) 

It  has  been  observed  by  Saxton  (1952a)  that  the  dielectric  behavior 
of  liquid  water  is  continuous  down  to  -8®C.  Hence  the  above  model  is 
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expected  to  be  equally  valid  for  liquid  water  at  temperatures  below 
0°C  (Ulaby,  1986). 

Figure  2  gives  the  static  dielectric  constant  (Eswo) 
relaxation  time  (Xsw)  aqueous  sodium  chloride  (NaCl)  solutions  as  a 

function  of  normality  and  temperature  as  computed  from  Eqs.  (19)  to 
(21),  It  is  interesting  to  note  in  Figure  2(a)  that  £5^-0  for  water  can 

drop  down  to  even  40  at  high  salinity  and  temperature.  It  can  be 

noted  from  Eq.  (5)  that  for  small  frequencies  (  in  the  MHz  range)  the 

conductivity  term  (i.e.,  the  second  term)  dominates  the  dielectric 
relaxation  term  (i.e.,  the  first  term),  but  at  higher  frequencies  (GHz 
range)  the  relaxation  effect  could  be  significant. 

The  models  described  earlier  in  this  section  are  applicable  to 
water  in  bulk.  However,  some  researchers  have  pointed  out  that  the 
dielectric  behavior  of  adsorbed  water  (i.e.,  water  in  contact  with  the 
surface  of  pores  of  rocks  and  soils)  is  slightly  different  from  that  of 
bulk  water  because  the  relaxation  frequency  for  adsorbed  water  is  a 
bit  lower  (Hoekstra,  1974)  and  the  ionic  conductivity  is  altered 

(Keller,  1966).  When  an  electrolyte  is  in  contact  with  a  surface, 

several  layers  of  water  molecules  will  become  adsorbed  to  the 
surface.  The  conductivity  of  water  in  this  oriented  adsorbed  phase  is 
higher  than  the  conductivity  of  free  water,  and  so  contributes  to  the 
overall  conductivity  of  a  rock  (Keller,  1966).  On  the  other  hand,  the 
increased  pressure  in  the  adsorbed  layers  increases  the  viscosity  of 
the  water  and  decreases  the  mobility  of  ions.  If  many  ions  are 
adsorbed,  the  water  conductivity  may  be  significantly  reduced 
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Temperature  (op)  Normality 

(a)  ^  (b) 

Figure  2.  (a)  Static  Dielectric  Constant  (b)  Relaxation  Time, 

of  Aqueous  NaCl  Solutions  (from  Stogryn,  1971) 


dUL K  WATER  RESiS  TIVITY  { ohm-m ) 

Figure  3.  Comparison  of  Resistivities  of  Pore  and  Bulk  Water 
(from  Keller,  1966) 
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(Keller,  1966).  Unfortunately,  a  mathematical  model  to  account  for 
the  above  effects  has  not  been  developed  so  far.  If  all  these  factors 
that  affect  the  dc  conductivity  of  water  in  pores  of  a  rock  are 

considered,  then  the  dc  conductivity  might  behave  as  shown  in 

Figure  3  (see  Keller,  1966).  The  linear  relationship  holds  if  there 

were  no  interaction  between  the  water  and  the  rock  framework.  The 
solid  curve  shows  the  more  probable  relation  between  rock 
resistivity  and  water  resistivity  when  interaction  between  the  water 

and  the  solid  minerals  takes  place.  At  high  salinities  (high 

conductivity  and  low  resistivity),  the  resistivity  of  the  saline  w'ater  in 
the  pores  is  increased  slightly  by  the  greater  viscosity  of  the  water 

adsorbed  on  the  grains.  At  lower  salinities  (lower  conductivities  and 
high  resistivities)  the  resistivity  of  pore  water  approaches  a 
maximum  value  rather  than  increasing  indefinitely.  This  limit  is 
determined  by  the  amount  of  interaction.  Both  the  increase  in 

resistivity  at  high  salinities  and  the  decrease  in  resistivity  at  low 
salinities  is  more  pronounced  in  fine-grained  rocks  than  in  coarse¬ 
grained  rocks  (Keller,  1966).  For  clay-rich  rock,  the  limiting  value 

may  be  as  low  as  0.1  ohm-m,  while  for  coarse  grained  rock  (e.g., 
concrete)  with  a  low  cation  exchange  capacity,  this  value  may  be  as 
high  as  10  ohm-m  (Keller,  1966),  which  is  the  value  shown  in  Figure 
3. 

The  salinity  of  pore  water  in  concrete  bridges  is  usually  much 

higher  than  10  ppt  (see  computations  in  the  chapter  on  'numerical 

studies').  For  this  lower  limit  of  S  =  10  ppt  and  T  =  20  ®C,  one 

s  w 

obtains  an  ionic  conductivity  of  1.5  mho/m  which  corresponds  to  a 
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resisiiviiy  of  G.67  ohm-n..  It  can  be  seen  from  Figure  3  that  li.is 
value  is  well  within  the  "linear"  range  for  a  coarse-grained  material. 
The  pore  water  resistivity  is  only  slightly  higher  in  this  range. 

The  above  discussion  holds  good  for  the  dc  case.  At  higher 
frequencies,  the  deviation  from  linear  relationship  is  further  reduced 
because  of  the  reduction  in  relaxation  frequency  for  adsorbed  water 
which  tends  to  increase  its  dielectric  conductivity  (see  Eqs.  5  and  13). 

2.3  Models  for  Dielectric  Properties  of  Concrete 

This  section  presents  three  models  which  can  be  used  to 
predict  the  dielectric  properties  of  a  mixture  from  the  properties  and 
volumetric  proportions  of  the  constituents.  These  models  have  been 
developed  and  verified  for  partially  saturated  rocks  and  will  now  be 
applied  to  concrete.  The  most  significant  constituents  of  concrete  are 
coarse  aggregate,  fine  aggregate,  cement  paste,  air,  water  and  salt. 
Reinforcing  bars  will  not  be  considered  here  as  a  constituent,  but 
their  effect  will  be  separately  superimposed  in  the  waveform 
synthesis  model  discussed  later.  It  has  already  been  mentioned  in 
the  previous  section  that  only  the  chlorides  which  are  actually 
present  in  solution  form  at  any  given  time  affects  the  dielectric 
properties  of  water  and  of  the  concrete  mixture  as  a  whole.  The 
coarse  aggregate  as  well  as  the  fine  aggregate  in  the  mortar  matrix 
have  a  relative  dielectric  constant  in  the  range  of  4.0  to  7.0  w-hich  is 
a  typical  value  for  common  aggregates.  The  porosity  of  the  coarse 
aggregate  used  in  New  England  area  (mainly  granite)  is  usually  very 
low  (about  1%),  hence  it  holds  very  little  water.  Thus,  the  aggregate's 
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dielectric  constant  has  a  very  small  imaginary  part  (loss  factor) 
which  may  be  ignored.  The  mortar  matrix  is  very  porous  and  the 
saline  water  in  these  pores  is  the  primary  constituent  which 
contributes  to  the  loss  factor  of  the  concrete  mixture.  Thus,  concrete 
can  be  visualized  as  a  three-phase  mixture  consisting  of  solid 
particles  and  air  which  have  a  real  dielectric  constant  and  saline 
water  which  has  a  complex  dielectric  permittivity.  Three  mixture 
models  which  are  applicable  to  concrete  are  presented  below. 


(a)  Complex  Refractive  Index  Method  ('CRIMl 
Refractive  index  (t]j)  for  a  medium  is  defined  as: 


Tjt  =  CoV^ltEoe 


(26) 


where, 

Cq  =  velocity  of  electromagnetic  waves  in  vacuum  = 

3  x  10^  m/sec 

Pt  =  permeability  of  the  medium  =  =  47cxl0’2  henry /m 


for  vacuum  and  most  dielectric  materials 

e  =  relative  complex  dielectric  permittivity  of  the  medium 

Cq  =  dielectric  permittivity  of  free  space  =  8.854  x  10*^2 
farad/m. 


The  complex  refractive  index  method  (CRIM)  asserts  that  the 
effective  complex-refractive  index  for  the  mixture  is  given  by  the 
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volume  average  of  the  complex  refractive  indices  of  the  constituents 
(Feng,  1985): 

Ve  =  (1-<|))  Ve^r+  (l-S)(j)Ver+  (27) 

where, 

<|>  =  porosity  of  concrete  =  (volume  of  voids)/(total  volume 

of  concrete) 

S  =  degree  of  saturation  =  (volume  of  water)/(volume  of 
voids) 

=  relative  dielectric  permittivity  of  concrete  solids  =  5.0 
(real) 

£a  ss  relative  dielectric  permittivity  of  air  h  1.0  (real) 

£sw  =  relative  complex  dielectric  permittivity  of  water 

(obtained  using  the  model  presented  in  the  previous 
section) 

£  =  relative  dielectric  permittivity  of  resulting  concrete 

mixture  (a  complex  quantity) 

The  CRIM  has  been  widely  used  for  soils,  rocks  and  concrete 
because  of  its  simplicity,  but  has  no  theoretical  basis.  It  gives 
reasonable  results  for  £"sw  -  1  but  gives  highly  inaccurate  results  for 

high  salinity  or  low  frequency;  that  is,  when  e"sw  »  1  (Feng,  1985). 
A  comparison  of  the  results  predicted  by  this  model  and  other  more 
sophisticated  models  for  the  range  of  properties  and  frequencies 
encountered  in  concrete  bridge  decks  has  been  carried  out  as  a  part 
of  this  research  and  is  presented  in  Chapter  4  of  this  report. 


W 
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Sometimes  CRIM  has  been  used  as  a  real  mixture  law  to  predict 
the  real  part  of  e  by  considering  only  the  real  part  of  dielectric 

permittivity  of  water  (e'sw)  and  ignoring  the  imaginary  part.  The 
accuracy  of  such  a  simplification  has  also  been  discussed  in  the 

chapter  on  numerical  studies  (Chapter  4). 

(b)  Continuous  Grain  Size  Distribution  Model 

This  model  assumes  that  the  solid  grains  and  air  molecules  in 

the  concrete  mixture  are  spherical  in  shape  and  have  a  continuous 

size  distribution.  The  model  has  been  derived  by  Feng  and  Sen 
(1985)  based  on  the  effective  medium  theory,  which  asserts  that  the 
effective  complex  relative  dielectric  constant  of  a  mixture  (e)  is  given 
by  (Feng,  1985): 


where  denotes  the  dielectric  permittivity  of  the  individual 
constituents  and  vj  the  corresponding  volume  fraction  of  the 

components.  It  has  already  been  mentioned  earlier  that  water  is  the 
only  conducting  phase  in  concrete,  and  the  solid  particles  and  air 
have  a  dielectric  permittivity  with  a  negligible  imaginary  part.  The 
above  effective  medium  theory  is  not  directly  applicable  to  concrete 
or  rocks  because  it  predicts  that  the  dc  conductivity  becomes  zero 
when  the  volume  fraction  <|»S  of  the  conducting  phase  (water)  drops 
below  a  critical  value  of  33%  (Feng,  1985).  For  most  rocks  and 
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concrete,  <|>S  is  almost  always  less  than  this  critical  value  (Feng, 
1985).  This  limitation  of  the  effective  medium  theory  can  be 
removed  by  following  a  systematic  procedure  in  which  one  starts 
with  a  background  matrix  of  water  in  which  small  amounts  of  solid 
particles  and  air  molecules  are  added  in  a  number  of  steps.  The 
smallest  particles  are  added  first  and  then  followed  by  larger 
particles.  Only  the  relative  sizes  of  these  particles  need  to  be 
considered  in  the  modelling.  After  each  small  addition  the  effective 
dielectric  permittivity  of  the  mixture  is  computed  using  Eq.  (28),  and 
the  resulting  mixture  serves  as  a  homogeneous  background  matrix 
for  the  next  phase.  The  additions  should  be  done  in  such  a  way  that 
the  final  mixture  has  the  constituents  in  the  desired  proportions. 
The  amount  of  addition  or  inclusion  of  particles  at  any  stage  should 
be  sufficiently  less  than  67%  of  the  volume  of  the  background  matrix 
from  the  preceding  stage.  This  removes  the  limitation  of  the 
equivalent  media  theory  (Feng,  1985). 

The  model  proposed  by  Feng  and  Sen  (1985)  assumes  a 
continuous  size  distribution  of  the  solid  particles  and  air  molecules, 
and  considers  infinitesimal  inclusion  of  these  panicles  at  each  step  as 
described  above.  It  can  be  shown  from  Eq.  (28)  that  a  change  de  in  e, 

the  complex  dielectric  constant  of  the  mixture  at  any  stage,  due  to 
inclusions  dv^  of  solid  grains  and  dv^  of  air  molecules  is  given  by 

(Feng,  1985): 

d£.  dVm  .  Cm  -  €  dVa  .  £a  >  £ 

3e  Vm  +  Vw  +  Va  Em  +  2£  +  Vw  +  Vj  Eg  +  2£ 


(29) 


-  25  - 


Here  v_,  v  and  v„  denote  the  volume  fraction  of  solid  particles, 
m  w  a 

water  and  air  in  the  mixture. 

In  general,  dv^  and  dv^  can  have  arbitrary  relation  at  each 

stage  of  the  mixing  process,  but  with  the  following  constraint  on  the 
final  ratio  of  the  constituents: 

''a  '  '"m  •  ''w  ~  ^  :  (1-  (|)):  S  0  (30) 


Feng  and  Sen  (1985)  assume  the  following  simple  relationship 
between  dv^^  and  dv^  at  each  step  of  the  mixing  process: 


^  (final)  =  I.' = 


dv 


m 


Vm 


(31) 


Thus,  the  medium  has  been  built  up  by  keeping  the  volume  of  water 
v^,  fixed  and  adding  infinitesimal  amounts  of  solid  grains  and  air 

molecules.  The  connectivity  of  water  molecules  is  ensured  at  each 
step  in  this  formulation  even  for  unsaturated  case,  which  is  very 
crucial  for  predicting  the  dielectric  conductivity.  Substituting  Eq. 
(31)  into  Eq.  (29)  and  integrating  the  resulting  equation  with  va/vm 
=  n'  held  fixed  and  using  the  boundary  condition  in  terms  of  the  loss 
factor  of  the  mixture  as  e"  =  e"sw  when  <|>S=1  (with  S  -»  I  being  taken 

first),  one  obtains  an  expression  for  the  effective  relative  dielectric 
constant  of  the  mixture  e  as  (Feng,  1985): 
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where, 

n’  =  4)  (l-S)/(l-4») 

a’  =  -2(1 +n') 

b’  =  2  e  -  e  +  n'  (2£„  -  e  ) 
m  a  ^  a 

c’  =  e.  (l-t-n') 
m  a 

m'  =  -[(2e_^  +  e^)  +  n'(2e^  +  E„)l  /  (2V  b'^  -  4a'c' ) 

In  the  dc  limit  (co  0),  Eq.  (32)  reduces  to  (Feng,  1985): 

=  Oj  S""  (t.™" ,  n"  =  m"  =  I  (33) 

where  and  refer  to  the  dc  conductivity  of  the  mixture  and 
water  respectively. 

Eq.  (33)  is  nothing  but  the  well  known  expression  for  Archie's 
Law  (McNeill,  1980)  which  is  commonly  used  to  predict  the  dc 
conductivity  for  fully  and  partially  saturated  soils  and  rocks.  A 
higher  value  of  n"  and  m"  (about  1.65  to  2.15)  has  been  reported  in 
the  literature  for  fully  saturated  rocks  (Keller,  1966).  A  value  of  2.0 
has  been  recommended  by  McNeil  (1980)  for  partially  saturated  soil 
based  on  experimental  results.  It  should  be  noted  that  Archie's  law 
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is  valid  only  in  the  dc  limit  whereas  Eq.  (32)  predicts  a  frequency 
dependent  behavior.  The  experimental  data  presented  by  Feng  and 
Sen  (1985)  shows  that  Eq.  (32)  gives  good  results  for  rocks  without 
flat  shaped  flaky  grains,  which  is  typical  of  claylike  materials.  The 
commonly  used  concrete  aggregates  do  not  contain  such  flaky  grains 
in  significant  amount,  and  hence  the  assumption  of  spherical  grain 
shape  is  a  reasonable  one. 

(c)  Discrete  Grain  Size  Distribution  Model 

The  previous  model  assumed  a  continuous  size  distribution  of 
spherical  grains.  On  the  other  hand,  a  model  proposed  by  Madden 
and  Williams  (1989)  using  the  equivalent  media  theory  considers  a 
discrete  particle  size  distribution  of  spherical  grains,  which  is  thus  a 
more  realistic  model  for  concrete.  This  model  is  also  built  in  steps 
with  the  smallest  particles  being  added  first  and  followed  by  larger 
particles.  As  mentioned  earlier,  only  the  relative  sizes  of  these 
particles  need  to  be  considered  in  the  modelling.  Grain  size 
distribution  is  an  important  factor  because  it  affects  the  internal 
structure  of  the  given  material,  which  in  turn  affects  its  dielectric 
properties.  Water  represents  the  smallest  particle  size  and  hence 
form  the  first  background  matrix,  and  an  iterative  procedure  is 
followed  to  recompute  the  dielectric  properties  after  each  addition  of 
a  small  amount  of  grains  to  this  matrix.  The  connectivity  of  water 
molecules  is  ensured  at  each  step  even  for  unsaturated  mixtures,  just 
as  in  the  previous  model. 


The  proposed  formulation  deals  with  complex  conductivities 
instead  of  complex  dielectric  constants,  but  these  two  concepts  can 
easily  be  interchanged.  Following  the  notation  of  Eqs.  (3)  to  (7),  the 
complex  relative  dielectric  constant  (e)  of  a  mixture  can  be  written  in 
terms  of  real  and  imaginary  parts  as: 


e  =  e’  +  ie 


•• 


or 


e  =  e’  + 


(OCo 


(34a) 

(34b) 


Here  a  denotes  the  dielectric  conductivity  of  the  medium  defined  as: 

o=:(oe"e^  (35  a) 

The  complex  conductivity  (o  )  is  defined  as: 

c 

=  a  -  i  e’£  0)  (35b) 

Co 

From  Eqs.  (34)  and  (35b),  the  relationship  between  £  and  a  can  be 

c 

obtained  as: 


a  =  -  i  0)  ££ 
c  o 


(36) 


One  can  also  formulate  this  model  using  the  complex  resistivity 
approach. 
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There  are  two  basic  approaches  to  the  equivalent  media 
calculations  (Madden,  1989).  One  is  the  symmetric  equivalent 
medium  theory,  which  assumes  all  particles  in  the  mixture  to  have 
spherical  grain  shape.  The  second  approach,  which  is  called  the 
unsymmetric  equivalent  medium  theory,  arises  from  the  important 
case  where  one  of  the  components  of  the  media  is  a  matrix  whose 
shape  is  essentially  unknown.  This  is  true  with  concrete  mixture 
where  the  water  phase  forms  the  background  matrix.  The 
previously  discussed  continuous  model  and  the  discrete  model  to  be 
discussed  in  this  section  both  fall  under  the  unsymmetric  equivalent 
medium  theory. 

Madden  and  Williams  (1989)  have  derived  the  following 
iterative  mixture  law  for  the  above  mentioned  discrete  case: 


On  =  On-i 


^  ^(On»l  ~  Oj)Pj  On 

^  (Gj  +  2an) 


(37) 


where, 

Gn  =  complex  conductivity  of  the  mixture  up  to  n^^^  step 

On-1  =  complex  conductivity  of  the  mixture  up  to  (n-l)^^  step 

Oj  =  complex  conductivity  of  inclusion 

Pj  =  volume  fraction  of  inclusion 
k'  =  number  of  inclusions  at  each  step. 


Eq.  (37)  is  a  quadratic  equation  in  On  which  can  be  solved  iteratively 
or  using  direct  closed  form  solution  at  each  step  to  obtain  two  roots 
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for  On-  If  there  is  only  one  inclusion  at  each  step  then  k’=l.  Among 
the  two  roots  of  Eq.  (37),  the  root  corresponding  to  the  positive  root 
of  the  discriminant  is  the  one  compatible  with  the  sign  convention 
used  in  Eq.  (35). 

The  inclusions  at  each  step  substitute  the  grains  in  the 
background  matrix  up  to  the  previous  step,  thus  altering  the 
proportion  of  these  grains  at  each  step.  Thus  the  volume  fractions 
(Pj)  should  be  back-calculated  in  such  a  way  that  the  matrix  in  the 

final  step  has  all  the  constituents  in  the  desired  proportion.  For  the 
case  of  one  inclusion  at  each  step,  the  following  recursive  relation  can 
be  used  to  compute  the  volume  fraction. 


Pi  = 


_ ^ _ 

(l-Pi)  (I-P2)  ...d-Pi-i) 


(38a) 


Pi  =Ai 


(38b) 


where, 

=  volume  fraction  of  the  inclusion  in  i^h  step 
Aj  5=  desired  final  proportion  of  the  inclusion  in  i^h  step 

^1  =  final  proportion  for  the  largest  size  particle. 

Among  the  three  mixture  models  presented  in  this  chapter, 
CRIM  is  the  simplest  one  and  has  no  theoretical  foundation.  The 
Continuous  model  and  the  Discrete  model  are  more  robust.  Since 
concrete  consists  of  discrete  particle  sizes,  the  Discrete  model  better 
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represents  the  physical  reality.  A  comparison  of  the  results 
predicted  by  the  three  mixture  models  will  be  presented  in  the  next 
chapter. 
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CHAPTER  3 

NUMERICAL  STUDIES 

3.1  Introduction 

This  chapter  presents  numerical  results  obtained  using  the 
models  and  equations  discussed  in  the  previous  chapter.  Firstly,  the 
data  used  for  this  numerical  study  and  its  sources  will  be  presented. 
The  subsequent  sections  will  describe  the  dielectric  properties  of 
concrete  obtained  by  means  of  the  various  models.  These  results  will 
then  be  examined  and  their  implications  will  be  discussed  at  the  end 
of  this  chapter. 

3.2  Data  Used  for  Numerical  Study 

The  various  input  parameters  in  these  models  for  dielectric 
properties  of  concrete  are  temperature,  frequency,  porosity  and 
degree  of  saturation  in  concrete,  salinity  of  the  pore  water  and 
dielectric  properties  of  concrete  solids.  In  case  of  the  discrete  grain 
size  distribution  model,  a  more  detailed  volumetric  proportions  of 
the  concrete  constituents  is  also  required. 

Temperatures  of  20  OC  and  5  have  been  considered  in  this 
study.  The  former  represents  conditions  during  fairly  warm  summer 
days  while  the  latter  is  more  typical  of  winter  months.  The  radar 
antennas  typically  used  for  concrete  bridge  decks  transmit  energy  in 
a  narrow  frequency  band  with  half  band-width  spread  of  about  0.6 
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to  1.4  GHz  and  a  central  frequency  of  1  GHz  (Maser  and  Halabe, 
1988).  The  porosity  of  newly  cast  concrete  is  usually  of  the  order  of 
0.1  (Popovics,  1985)  but  disintegration  due  to  freeze-thaw  cycles  can 
cause  concrete  in  actual  bridge  decks  to  have  a  higher  porosity  of 
about  0.15,  and  even  upto  0.2  in  some  extreme  cases  (Cady,  1983). 
For  the  purpose  of  this  numerical  study,  porosities  of  0.10  and  0.15 
have  been  considered.  The  degree  of  saturation  of  the  top  and 
bottom  one  inch  layer  of  concrete  in  conventional  concrete  bridge 
decks  and  the  top  one  inch  of  stay-in-place  (SIP)  concrete  decks 
(which  are  covered  from  bottom)  varies  between  40  or  50%  in  warm 
summer  months  to  75%  in  cold  winter  months  (Carrier,  1975).  The 
inner  portion  of  the  concrete  usually  retains  a  higher  degree  of 
saturation,  fluctuating  between  70%  and  95%  round  the  year  (Carrier, 
1975).  This  numerical  study  presents  results  for  the  entire  range  of 
0  to  100%,  but  more  emphasis  has  been  placed  on  the  50%  to  100% 
range. 

No  direct  measurement  of  salinity  of  water  present  in  the  pores 
of  concrete  exposed  to  natural  environment  has  been  found  in  the 
literature.  One  exception  to  this  is  the  case  of  concrete  exposed  to 
marine  environment,  which  is  not  the  area  of  application  considered 
here.  Hence,  the  salinity  of  pore  water  in  concrete  bridge  decks  has 
been  indirectly  ascertained  based  on  field  measurements  of  other 
related  parameters.  Values  for  free  chloride  and  moisture  content 
(expressed  as  percentage  by  weight  of  concrete)  in  core  samples 
taken  from  deteriorated  bridges  in  New  England  areas  are  presented 
in  a  report  by  Maser  (1989).  In  two  cases  typical  of  older  (>  25 
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years)  decks,  the  free  chloride  content  in  top  one-inch  layer  of 
concrete  was  found  to  be  0.25%  and  0.34%  (by  weight  of  concrete) 
and  the  corresponding  moisture  content  was  4.8%  and  3.8%  (by 
weight  of  concrete).  To  convert  the  percentage  chloride  content  to 
NaCl  equivalent,  a  multiplication  factor  equal  to  the  ratio  of 
molecular  weight  of  NaCl  (=  23.0  +  35.5  =  58.5)  to  the  atomic  weight 
of  chlorine  (=  35.5)  has  to  be  applied.  This  ratio  turns  out  to  be  1.65. 
Dividing  the  NaCl  equivalent  of  the  free  chloride  content  by  the 
corresponding  moisture  content,  the  salinity  of  pore  water  in  the  two 
cases  is  obtained  as  86  and  148  parts  per  thousand  by  weight  (ppt). 
A  similar  calculation  on  another  deteriorated  bridge  deck  showed  a 
salinity  of  132  ppt  in  the  top  one  inch  layer  of  concrete.  The  chloride 
content  reduces  significantly  with  depth  in  the  concrete  deck  (Maser, 
1989;  Cady,  1983).  Data  presented  by  Cady  and  Weyers  (1983) 
shows  that  the  salinity  is  usually  21  ppt  for  new  concrete  decks. 
Such  a  small  amount  of  salinity  is  inherent  in  concrete  because  of 
impurities  in  mixing  water.  Sometimes,  salt  is  added  to  concrete 
during  the  mixing  process  to  reduce  the  setting  time.  For  the 
purpose  of  this  numerical  study,  salinity  values  ranging  between  12 
and  90  ppt  have  been  considered.  This  range  represents  the  salinity 
of  a  typical  deck  during  its  first  20  years  of  service  life. 

The  dielectric  properties  of  saline  water  can  be  computed  using 
the  equations  presented  in  Section  2.2.  Now,  the  only  remaining 
parameters  are  the  dielectric  constant  for  air,  which  is  1.0  and  the 
dielectric  constant  for  concrete  solids,  which  is  taken  as  5.0.  Typical 
values  for  granite,  sandstones  and  limestones  usually  range  between 
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4.0  and  7.0.  The  imaginary  pan  of  the  complex  dielectric 
permittivity  is  zero  for  air,  and  negligible  for  concrete  solids- 
compared  to  the  value  for  saline  water. 

The  discrete  grain  size  distribution  model  requires  more 
detailed  volumetric  proportions  of  the  concrete  constituents.  The 
amount  of  water  in  the  mix  is  given  by  <t)S,  and  the  remaining  void 
volume  equal  to  (1-S)<|>  is  filled  with  air.  About  two-thirds  of  this  air 
volume  is  in  the  form  of  fine  pores  of  size  less  than  100  nm  and  the 
remaining  one-third  constitutes  coarse  pores  of  approximate  size 
0.05  mm  (Neville,  1983).  The  fine  pores  are  the  result  of  moisture 
evaporation  during  the  setting  process,  while  the  coarse  pores  are 
entrained  ai>  generated  during  the  mixing  process.  The  entrained  air 
helps  in  reducing  freeze-thaw  damage  of  exposed  concrete  (Neville, 
1983).  The  fine  aggregates  (i.e.,  mortar  paste  of  cement  and  sand) 
usually  consists  of  panicles  with  size  less  than  0.125  mm  (Neville, 
1983).  The  coarse  aggregates  (gravels)  in  a  typical  concrete  mix 
usually  have  a  size  range  of  1/2  to  3/2  inch  (Neville,  1983)  and  a 
volumetric  content  of  about  50%.  The  remaining  half  of  the  volume 
consists  of  fine  aggregates  and  pores  filled  with  air  and  water. 

As  mentioned  earlier,  the  discrete  grain  size  model  is  built  in 
steps  with  the  smallest  particles  being  added  first  and  followed  by 
larger  particles.  At  each  step  the  volume  fraction  of  the  inclusion 
(Pj)  should  be  significantly  less  than  67%  of  the  total  volume  of 

background  matrix.  To  take  care  of  this  constraint,  the  volume  of 
particles  in  each  category  was  further  subdivided  into  small  bins  as 
follows: 


-  36  - 


(a)  volume  fraction  of  water  =  (j»S  =  one  bin  (Aji) 

2  _ 

(b)  volume  fraction  of  air  in  fine  pores  =  ^  (1'S)<1)  =  3  bins  (Ag, 
A9,  Aio) 

(c)  volume  fraction  of  fine  aggregates  =  (0.5  -  0)  =  3  bins  (A5, 
Ae,  A7) 

(d)  volume  fraction  of  air  in  coarse  pore  =  ^(1-S)(|>  =  1  bin  (A4) 

(e)  volume  fraction  of  coarse  aggregates  =  0.5  =3  bins  (Aj,  A2, 
A3) 

The  final  volume  fraction  in  the  above  bins  have  been  chosen  as 
follows: 

All  = 

Aio  =  i  (1  -  S)it. 

A9  =  i  (l-S)* 

As  =  i  (1-5)0 
At  =  i  (0.5  -  0) 

As  =  i  (0.5  -  $) 

As  =  ^  (0.5  -  d)) 
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A4  =  ^  (1-8)4) 

A3  =  A2  =  Ai  = 

3 

While  building  the  mixture  in  this  numerical  model,  the  smallest  size 
particle  (i.e.,  Ajj)  is  added  first  and  the  other  particles  Ajq  Aj 

are  added  sequentially  with  one  bin  at  each  step.  It  may  be 
mentioned  here  that  the  final  complex  dielectric  permittivity  of  the 
concrete  mixture  predicted  by  this  model  is  not  very  sensitive  to  the 
assumptions  made  above  regarding  the  volumetric  content  of  the 
concrete  constituents  and  the  volume  distribution  in  the  bins,  as  long 
as  all  the  Pj’s  are  significantly  less  than  0.67,  as  will  be  shown  in  the 

next  section. 

3.3  Numerical  Results 

This  section  presents  the  numerical  results  obtained  using  the 
models  and  equations  presented  in  Chapters  2  and  3  and  the  data 
given  in  the  previous  section.  The  implications  of  these  results  are 
discussed  in  the  following  paragraph. 

(a)  Comparison  of  Complex  CRIM  vs  Real  CRIM  Model 

Table  1  shows  the  comparison  of  the  results  for  real  part  of 
concrete  dielectric  permittivity  predicted  using  CRIM  (Eq.  27)  and  a 
simplified  version  of  CRIM  which  has  been  termed  here  as  "Real 
CRIM".  The  CRIM  model  in  Eq.  (27)  uses  the  complex  dielectric 
permittivity  of  saline  water  (Ejw)  and  predicts  complex  permittivity 
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TABLE  1  COMPARISON  OF  COMPLEX  CRIM  vs  REAL  CRIM 

T  =  20OC,f=  1  GHz,<j>  =  0.1. 
e'(air)  =  1.0,  e*  (solids)  =  5.0 


Degree 

of 

1  ^sw  “  PP^*  ^  sw  ~  77.0 

Sgw  ~  ^9  PP^' 

£  sw  ~  57.4 

Saturation 

S 

e’  Using 
CRIM 

e'  Using 

Real  CRIM 

e'  Using 
CRIM 

e'  Using 

Real  CRIM 

0.0 

] 

*  4.5 

1 

4.5 

4.5 

4.5 

0.2  1 

1  5.2 

5.1 

5.4 

5.0 

0.4  j 

] 

i  5.9 

5.9 

6.3 

5.6 

0.6 

] 

6.7 

6.7 

7.2 

6.3 

0.8  1 

j 

1  7.6 

j 

7.5 

8.2 

7.0 

1.0 

■  8.5 

8.4 

9.2 

7.7 
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(e)  for  concrete.  The  Real  CRIM  model  uses  the  real  part  of  water 
permittivity  (e’sw)  predicts  the  real  part  of  concrete  permittivity 
(e').  This  real  version  of  CRIM  is  commonly  used  because  of  its 
simplicity.  It  can  be  seen  from  Table  1  that  these  two  versions  of 
CRIM  model  predict  almost  the  same  results  for  low  salinity  of  water 
(Ssw  =  12.0  ppt)  but  deviate  to  some  extent  for  high  salinity  (Sgw  = 
80.0  ppt).  The  phase  velocity  and  the  reflection  and  transmission 
coefficients  involve  the  square  root  of  real  dielectric  permittivity 
(V?),  and  the  Real  CRIM  predicts  a  value  lower  by  about  9%  for  V?  in 
the  high  salinity  and  degree  of  saturation  range,  as  can  be  seen  from 
Table  1.  Such  a  discrepancy  is  small  in  making  velocity  calculations, 
but  can  be  significant  in  comparing  reflection  coefficients.  Also,  the 
Real  CRIM  model  does  not  predict  the  loss  factor  which  is  related  to 
attenuation. 

(b)  Comparison  of  the  Three  Mixture  Models  for  Concrete 

Table  2  presents  the  results  predicted  by  the  three  complex 
mixture  models  for  various  values  of  porosity,  salinity  and  degree  of 
saturation  at  the  1  GHz  frequency.  The  Discrete  model  could  not 
predict  any  value  for  S  =  0.0  because  for  this  case  the  volume 
fraction  of  the  inclusions  (Pj)  happen  to  exceed  the  critical  value  of 

67%  due  to  absence  of  water.  It  is  possible  to  avoid  this  problem  by 
changing  the  bin  sizes  for  air  when  S  =  0.0,  which  was  not  done  here. 
It  can  be  seen  from  this  table  that  the  real  part  of  concrete  dielectric 
permittivity  predicted  by  the  three  models  are  very  close  (within  5% 
for  for  all  cases.  Field  measurements  on  actual  bridge  decks 
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TABLE  2  COMPARISON  OF  THE  THREE  MIXTURE  MODELS 
(T  =  20  OC,  f  =  1  GHz,  £’(air)  =  1.0,  e*  (solids)  =  5.0) 


(a)  (|) 

—  0.1,  Ssw  - 

12  ppt 

Degree 

of 

Saturation 

CRIM 

Model 

Continuous 

Model 

Discrete 

Model 

S 

e’ 

ki(m-l) 

c* 

ki(m-l) 

e’  ki(m'l) 

0.0 

4.5 

0 

4.5 

0 

i 

1 

i  . 

0.2 

] 

5.2 

0.9 

5.2 

1.0 

i  5.2  0.9 

0.4 

5.9 

1.7 

6.0 

2.0 

5.9  1.6 

0.6 

6.7 

2.6 

6.8 

3.0 

6.6  2.3 

0.8 

7.6 

3.4 

7.6 

4.0 

7.3  3.0 

1.0 

8.5 

4.3 

8.5 

5.0 

7.9  3.7 

TABLE  2  (continued) 


(b)  <^  = 

O'L  Sjvv  - 

80  ppt 

Degree 

of 
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CRIM 

Model 

Continuous 

Model 

Discrete 

Model 

S 

e’ 

ki(m-l) 

1 

e* 
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E' 
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4.5 

0.0 

4.5 

0.0 

0.2 

5.4 

3.5 
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3.4 
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1 

6.2 

1 
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7.2 
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6.9 

13.0 

6.9 

9.3 

0.8  1 
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14.1 

7.6 

18.0 

7.6 

12.4 

1 

1.0 

9.2 

17.6 

8.4 

23.0 

8.3 

15.7 
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TABLE  2  (continued) 
(c)  (j)  =  0.1,  Sjvv  =  52  ppt 


Degree 

of 

Saturation 

CRIM 

Model 

Continuous 

Model 

Discrete 

Model 

S 

e’  ki(m'l) 

i 

e‘ 

ki(m-l) 

£’ 

ki(ni*l) 
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0.6 

7.0 
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0.7 

7.4 

9,5  I 

7.2 
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7.1 
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0.8 
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10.8 

7.6 
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7.5 

9.4 

0.9  1 

8.4 

12.2 

8.0 

15.0 

7.8 

10.6 

1.0 

8.8 

13.5 

8.4 

17.0 

8.1 

11.9 
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TABLE  2  (Continued) 
(d)  <}>  =  0.15,  Ssw  =  52  ppt 


Degree 

of 

Saturation 

CRIM 

Model 

Continuous 

Model 

Discrete 

Model  ' 

1 

S 

e'  ki(m*l) 

e* 
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E‘ 
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12.2 

7.7 
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8.3 

14.3 
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9.4 
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16.4 

0.9 
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9.6 

23.0 

9.4 

18.5 

1.0 
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20.3 

10.2 

25.0 

10.0  20.5 
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have  indicated  that  dielectric  permittivity  is  about  9.0  for 
deteriorated  concrete,  which  suggests  that  the  porosity  should  be 
around  0.15  (see  Table  2d).  There  is  a  significant  discrepancy  in  the 
value  of  attenuation  coefficient  (kj)  predicted  by  the  three  models. 

The  discrepancy  between  the  Discrete  model  and  the  CRIM  model  is 
about  13%  for  high  salinity  and  degree  of  saturation  with  the  CRIM 
Model  predicting  a  higher  attenuation.  The  Continuous  model 
predicts  an  attenuation  coefficient  which  is  higher  by  even  25  to  40% 
compared  to  that  predicted  by  the  Discrete  model  for  high  values  of 
salinity  (Sg^v  =  80.0  ppt)  and  degree  of  saturation.  Thus,  it  can  be 

observed  that  the  CRIM  model  is  closer  to  the  Discrete  model  than  to 
the  Continuous  model. 

The  effect  of  the  discrepancie."  in  the  attenuation  coefficient 
predicted  by  the  three  models  on  the  predicted  waveforms  can  be 
understood  from  the  following  bridge  deck  calculations.  Since  the 
Continuous  and  Discrete  models  predict  the  highest  and  lowest  values 
for  attenuation  coefficient  (kj),  we  will  consider  these  two  models 

here.  For  the  purpose  of  this  calculation,  let  us  consider  a  6  inch 
thick  concrete  deck  with  the  top  one  inch  layer  having  a  salinity  of 
80  ppt  and  degree  of  saturation  of  0.6,  and  the  bottom  five  inches  of 
concrete  with  a  salinity  of  12  ppt  and  degree  of  saturation  of  0.8. 
From  Table  2,  the  corresponding  attenuation  coefficients  using  the 
Continuous  model  are  13.0m"l  for  top  one  inch  layer  and  4.0m‘^  for 
the  bottom  five  inch  layer.  The  Discrete  model  predicts  values  for 
the  attenuation  coefficients  of  9.3m*  1  and  3,0m*  1,  respectively. 
Since  the  transmitted  wave  has  to  travel  twice  the  layer  thickness 
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before  it  is  received,  the  radar  vt^ave  amplitude  attenuates  by  a 
factor  of  exp(-2kix)  as  cau  be  seen  using  Eq.  (44),  Hence,  the 

attenuation  factor  of  wave  reflected  from  interface  between  the  two 
concrete  layers  is  obtained  as  0.52  using  the  Continuous  model  and 
0.62  using  the  Discrete  model.  The  attenuation  factor  for  the  bottom 
reflection  is  a  product  of  the  attenuation  factors  in  the  top  one  inch 
and  bottom  five  inches  of  concrete,  and  is  obtained  as  0.19  for 
Continuous  model  and  0.29  for  Discrete  model.  Thus,  the  discrepancy 
in  the  two  models  ranges  from  20%  to  50%.  If  one  includes  all  the 
reflection  and  transmission  coefficient  at  the  various  interfaces  (e.g., 
air-asphalt  and  asphalt-concrete  boundary)  then  the  discrepancy  in 
the  final  synthesized  waveform  using  the  two  models  could  be  either 
higher  or  lower  depending  upon  how  close  the  properties  of  adjacent 
layers  are.  If  one  assumes  that  the  attenuation  coefficient  (kj)  for 

the  bottom  five  inches  is  also  as  high  as  that  for  the  top  one-inch 
layer  (which  could  be  the  case  for  a  severely  deteriorated  concrete 
deck),  then  the  Continuous  and  Discrete  models  predict  an 
attenuation  factor  for  bottom  reflection  of  0.02  and  0.06, 
respectively. 

These  predicted  attenuations  are  consistent  with  field 
observations.  The  field  data  from  some  deteriorated  concrete  bridge 
decks  in  New  England  area  shows  that  under  certain  conditions  it  is 
difficult  to  see  reflection  from  the  bottom  of  the  decks  (Maser,  1989). 
The  bottom  reflections  predicted  here  are  more  representative  of 
field  data  than  previous  predictions  (see  Maser  and  Halabe,  1988). 
This  was  one  of  the  motivations  for  this  study  (see  Chapter  1). 
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However,  the  predictions  from  the  two  models  (Discrete  and 
Continuous)  are  significantly  different.  Of  the  three  models,  the 
Discrete  model  represents  the  closest  representation  of  reality  for 
reinforced  concrete.  The  Continuous  model  would  be  more 
appropriate  for  materials  with  a  continuous  size  distribution.  The 
Continuous  model  appears  to  predict  higher  conductivity  than  the 
Discrete  model.  This  can  be  understood  physically,  as  discussed  in 
the  next  paragraph.  The  CRIM  model,  which  has  no  physical  basis, 
predicts  values  somewhere  in  between  those  of  the  two  physical 
models. 

The  above  numerical  results  correspond  to  a  frequency  of  1 
GHz.  For  low  values  of  frequency  (MHz  range)  or  extremely  high 
values  of  salinity  the  loss  tangent  for  concrete  (tan  6  =  e’Ve’)  is  much 
greater  than  unity  (see  Eq.  46)  and  the  CRIM  model  which  has  no 
theoretical  basis,  breaks  down.  On  the  other  hand,  the  Continuous 
and  Discrete  models  are  valid  even  for  tan  5  much  larger  than  unity. 
Calculations  carried  out  under  this  study  (but  not  reported  here) 
show  that  the  dielectric  conductivity  of  concrete  (defined  by  Eq.  35a) 
is  usually  independent  of  frequency  for  frequencies  below  the  MHz 
range.  This  value  represents  the  dc  conductivity  of  the  medium,  and 
can  be  obtained  using  the  Archie’s  law  expression  for  Continuous 
model  with  an  exponent  of  1.5  (see  Eq.  33).  Additional  calculations 
have  shown  that  a  similar  equation  holds  for  the  Discrete  model  as 
well,  but  with  an  exponent  of  about  1.7.  The  product  <j)S  in  Eq.  (33) 
always  has  a  value  less  than  unity.  Thus,  the  Discrete  model  predicts 
lower  values  for  dc  and  dielectric  (high  frequency)  conductivity  and 
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hence  lower  attenuation  as  compared  to  the  Continuous  model.  This 
is  due  to  the  increased  connectivity  of  water  molecules  when  the 
grain  size  distribution  is  continuous.  If  the  number  of  bins  in  the 
Discrete  model  is  increased,  then  the  results  converge  to  the  values 
predicted  by  the  continuous  model. 

The  discrepancy  in  the  attenuation  coefficient  predicted  by  the 
three  models  must  be  considered  in  the  context  of  the  possible  error 
in  the  estimation  of  salinity  of  the  pore  water  in  concrete.  It  has 
already  been  mentioned  earlier  that  this  salinity  is  indirectly 
estimated  by  dividing  the  free  chloride  content  in  the  concrete 
(obtained  from  chloride  content  measurements  in  core  samples)  by 
the  moisture  content  of  the  concrete  (obtained  from  moisture  content 
determination  in  the  core  sample).  However,  it  is  unknown  whether 
the  salinity  of  the  water  increases  as  the  concrete  loses  its  moisture 
or  remains  the  same  with  change  in  degree  of  saturation.  No 
research  data  is  available  on  the  amount  of  free  chlorides  actually 
remaining  in  dissolved  state  in  concrete  bridge  decks.  A  numerical 
study  was  conducted  as  a  part  of  this  research  using  the  mixture 
models  described  earlier  in  which  the  salinity  of  water  was 
proportionately  increased  with  reduction  in  degree  of  saturation, 
thus  assuming  that  all  the  chlorides  remained  in  dissolved  state.  For 
a  degree  of  saturation  of  0.7,  this  increased  salinity  corresponds  to 
about  25%  increase  in  the  attenuation  coefficient  over  the  case  for 
which  the  salinity  is  kept  at  the  same  level  as  that  for  100% 
saturation,  regardless  of  the  model  used.  The  change  in  salinity  has  a 
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very  small  effect  on  the  real  pan  of  concrete  permittivity  (e’)  which 
can  also  be  seen  from  Table  2. 

In  the  numerical  study  of  the  Discrete  model,  a  particular  size 
distribution  has  been  assumed  for  the  concrete  constituents.  A 
numerical  study  was  conducted  to  determine  the  sensitivity  of  the 
predicted  results  to  the  size  distribution.  A  modified  concrete  with  a 
slightly  different  size  distribution  was  assumed  for  this  purpose.  The 
modified  concrete  was  assumed  to  have  65%  coarse  aggregates,  with 
the  remaining  35%  of  the  volume  occupied  by  fine  aggregates,  air 
pore  and  water,  as  opposed  to  the  50-50  ratio  in  case  of  the  original 
concrete.  In  the  original  concrete,  two-thirds  of  the  total  air  volume 
was  assumed  to  be  in  fine  pores  and  one-third  in  coarse  pores,  but  in 
the  modified  concrete,  the  total  air  volume  was  equally  distributed 
among  the  coarse  and  fine  pores.  The  porosity,  water  salinity  and 
degree  of  saturation  was  kept  at  the  same  level  for  both  cases.  Table 
3  shows  the  results  predicted  by  the  Discrete  model  for  the  original 
and  modified  concrete.  As  explained  earlier,  no  values  could  be 
obtained  for  S  =  0.0,  because  due  to  the  absence  of  water,  the  volume 
fraction  of  subsequent  inclusions  (Pj)  exceed  the  critical  value  of 
0.67. 

The  results  presented  in  Table  3  show  the  differences  in 
dielectric  properties  between  the  original  and  modified  concrete. 
The  real  part  of  dielectric  permittivity  (e’)  is  almost  the  same  in  the 
two  cases,  but  the  attenuation  coefficient  differ  by  as  much  as  8%. 
Thus,  it  can  be  concluded  that  the  attenuation  coefficient  varies 
significantly  with  the  size  distribution. 
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TABLE  3  EFFECT  OF  SLIGHT  CHANGE  IN  CONCRETE 
COMPOSITION  ON  RESULTS  PREDICTED 
BY  DISCRETE  MODEL 

T  =  20OC.  f  =  1  GHz,  <j)  =  0.1 
e'(air)  =  1.0,  e‘  (solids)  =  5.0 
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Since  the  Discrete  model  better  represents  the  true  concrete 
composition  than  the  Continuous  model,  it  should  lead  to  better 
predictions  for  the  electromagnetic  properties  of  concrete.  The 
theoretical  models  presented  here  are  yet  to  be  verified  by  means  of 
either  radar  field  data  from  actual  bridge  decks  or  small  scale 
laboratory  experiments. 

(c)  Study  of  Dielectric  Properties  of  Concrete  Using  Discrete 
Model 

This  section  investigates  the  typical  behavior  of  concrete 
dielectric  properties  as  a  function  of  the  various  input  parameters.  It 
has  already  been  mentioned  earlier  (see  Table  2)  that  the  salinity  of 
pore  water  has  a  very  small  effect  on  the  real  part  of  concrete 
permittivity  but  has  a  large  effect  on  the  loss  factor  and  attenuation 
coefficient.  The  frequency  dependence  of  dielectric  properties  of 
concrete  has  been  studied  using  the  Discrete  model,  and  the  results 
are  shown  in  Figures  4,  5,  6  and  7.  Figure  4  shows  that  the  real  part 
of  complex  concrete  permittivity  (e’)  is  not  very  sensitive  to 
frequency  in  the  range  of  0.6  to  3.0  GHz.  Thus,  the  wave  velocity 
(which  is  inversely  proportional  to  V?)  is  non-dispersive  in  this 
frequency  range.  Figures  5,  6  and  7  show  the  loss  factor,  attenuation 
coefficient  and  dielectric  conductivity,  and  reveal  a  stronger 
dependence  on  frequency.  The  dispersive  characteristic  is  stronger 
for  higher  degree  of  saturation,  because  water  is  the  component 
which  makes  concrete  a  dispersive  medium.  It  can  also  be  seen  that 
the  attenuation  coefficient  and  dielectric  conductivity  at  high 
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Figure  7.  Dielectric  Conductivity  of  Concrete  (Discrete  Model ; 
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T  =  20°C,  f  =  1  GHz,  e'  (air)  =  1.0,  e'  (solids)  =  5.0 
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Figure  8.  Permittivity  of  Concrete  as  a  Function  of  Degree 
of  Saturation  (Discrete  M ocel ) 
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frequency  (near  1  GHz)  is  much  higher  than  in  the  dc  limit  case.  An 
interesting  feature  that  can  be  noted  from  Figures  6  and  7  is  that  the 
frequency  dependence  of  the  attenuation  coefficient  and  the 
dielectric  conductivity  is  almost  linear  for  frequencies  above  0.6  GHz. 

The  discrete  model  used  here  can  not  compute  results  for  S=0 
due  to  numerical  problems  outlined  earlier.  Hence,  the  values 

corresponding  to  S=0  in  Figures  4,  5,  6  and  7  have  been  computed 
using  GRIM.  Note  that  GRIM  and  Gontinuous  model  both  predict  the 
same  results  at  S=0.  Figures  8  and  9  correspond  to  data  presented  in 
Table  2,  and  show  that  the  real  part  of  concrete  permittivity  (s')  and 
the  attenuation  coefficient  are  linearly  related  to  the  degree  of 
saturation.  This  suggests  that  the  complex  models  presented  in  this 
research  could  be  approximated  by  simpler  equations  for  the 

prediction  of  concrete  dielectric  properties  in  the  range  of  interest. 
This  suggestion  will  be  investigated  in  the  next  phase  of  this 
research. 

All  the  above  numerical  study  was  conducted  for  a 

temperature  of  20  ®G.  It  has  already  been  mentioned  earlier  (see 
Figures  1  and  2)  that  the  dielectric  properties  of  water  have  a 

significant  temperature  dependence.  A  numerical  study  was 
conducted  to  understand  the  effect  of  temperature  on  the  dielectric 
properties  of  concrete  and  the  results  have  been  tabulated  in  Table 
4.  It  can  be  seen  that  the  temperature  change  has  negligible  effect 
on  the  real  part  of  concrete  dielectric  permittivity  (e').  On  the  other 
hand,  the  attenuation  coefficient  reduces  by  about  25%  when  the 
temperature  drops  from  20  ©G  to  5  OG.  This  is  because  the  water 
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TABLE  4  EFFECT  OF  TEMPERATURE  CHANGE  ON 
DIELECTRIC  PROPERTIES  OF  CONCRETE 

Ssw  =  52.0  ppt,  f  =  1  GHz,  (|)  =  0.15 
e'  (air)  =  1.0,  £*  (solids)  =  5.0 


Degree 

of 

Saturation 

S 

T  =  20  OC 

T 

=  5  OC 

e’ 

ki(m-^) 

£’ 

ki(m-l) 

1 

1 

0.0  1 

1 

4.2 

0.0 

4.2 

0.0 

0.2 

5.4 

4.1 

5.3 

3.2 

0.4 

6.6 

8.2 

6.5 

6.3 

0.6 

7.7 

12.3 

7.6 

9.3 

0.8 

8.9 

16.4 

8.8 

12.4 

1.0 

10.0 

20.5 

10.0 

15.5 
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dielectric  conductivity  drops  significantly  with  temperature.  Field 
data  from  some  deteriorated  bridge  decks  collected  in  wa'm  summer 
days  have  indeed  shown  significantly  higher  attenuation  compared 
to  the  data  collected  from  the  same  bridge  decks  in  cold  winter  days 
(Maser,  1989). 


I 

I 
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CHAPTER  4 

CONCLUSIONS 

The  objective  of  this  research  was  to  characterize  the  velocity 
and  attenuation  of  electromagnetic  waves  in  reinforced  concrete  as  a 
function  of  frequency,  temperature,  moisture  content,  chloride 
content  and  concrete  mix  constituents.  This  objective  has  been  met 
by  adaptation  of  theoretical  models  which  have  been  developed  and 

verified  for  partially  saturated  rocks.  These  models  predict  the 
electromagnetic  properties  of  concrete  by  aggregating  the 
electromagnetic  properties  of  its  most  significant  constituents  which 
are  water,  salt,  aggregates  and  air.  Although  the  volume  fraction  of 
water  in  a  concrete  mixture  is  small,  it  has  the  most  significant  effect 
on  concrete  dielectric  properties  because  of  its  high  complex 

dielectric  permittivity.  The  presence  of  salt  in  the  water  further 
influences  the  concrete  dielectric  properties. 

Three  models  have  been  investigated  in  this  study.  One  is  a 
simple  three-phase  mixture  law  called  CRIM  which  has  no  theoretical 
foundation.  The  other  two  models  are  more  robust  and  have  a 
theoretical  basis.  The  Continuous  grain  size  distribution  model 
assumes  that  the  solid  grains  and  air  molecules  in  the  concrete 

mixture  are  spherical  in  shape  and  have  a  continuous  size 

distribution.  On  the  other  hand,  the  Discrete  grain  size  distribution 
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model  considers  a  discrete  particle  size  distribution  of  spherical 
grains,  which  is  thus  a  more  realistic  model  for  concrete. 

A  comparison  of  the  results  predicted  by  the  Discrete  model 
with  the  complex  CRIM  model  (Eq.  27)  shows  that  the  two  models 
agree  very  well  for  real  part  of  concrete  dielectric  permittivity  (e'). 
The  attenuation  coefficient  predicted  by  the  complex  CRIM  model  is 
slightly  higher  (within  15%)  compared  to  the  Discrete  model.  The 
Continuous  model  predicts  almost  the  same  value  for  e’  but  a 
significantly  higher  value  for  attenuation  coefficient.  This  is  due  to 
the  increased  connectivity  of  water  molecules  when  the  grain  size 
distribution  is  continuous. 

Numerical  studies  have  been  carried  out  using  the  Discrete 
model.  This  model  was  chosen  as  the  most  accurate  representation 
of  real  concrete.  The  results  have  shown  that  the  salinity  of  pore 
water  has  a  negligible  effect  on  the  real  part  of  complex  concrete 
permittivity  but  has  a  large  effect  on  the  loss  factor  and  attenuation 
coefficient.  The  real  part  of  complex  concrete  permittivity  (e’)  is  not 
very  sensitive  to  frequency  in  the  range  of  0.6  to  3.0  GHz.  Thus,  the 
wave  velocity  (which  is  inversely  proportional  to  V?)  is  non- 
dispersive  in  the  high  frequency  range.  However,  the  loss  factor, 
attenuation  coefficient  and  dielectric  conductivity  show  a  stronger 
dependence  on  frequency.  The  frequency  dependence  of  attenuation 
coefficient  and  dielectric  conductivity  is  almost  linear  for  frequencies 
in  the  0.6  to  3.0  GHz  range.  The  real  part  of  concrete  permittivity 
and  the  attenuation  coefficient  also  show  a  linear  dependence  with 
respect  to  the  degree  of  saturation  of  the  concrete  mixture.  This 
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suggests  that  the  complex  models  presented  in  this  study  can  be 
approximated  in  the  range  of  interest  by  simple  equations  which  can 
be  used  with  a  pocket  calculator.  Preliminary  numerical  studies 
have  shown  that  even  a  real  CRIM  model  (a  simpler  version  of  the 
complex  CRIM)  with  proper  modification  of  the  exponent  can 
approximate  the  other  two  complex  models  for  the  prediction  of  £' 
within  5  to  10%  in  the  high  frequency  range.  The  complex  part  of 
concrete  permittivity  (e")  which  is  related  to  the  dielectric 

conductivity  (o)  can  also  be  predicted  within  an  error  of  10% 
compared  to  the  Discrete  model  by  suitably  modifying  the  exponents 
in  the  simple  Archie's  law  (Eq.  33)  or  the  complex  CRIM  (Eq.  27). 
Hence,  future  research  should  focus  on  developing  these  simpler 

equations  to  approximate  the  complicated  models  for  the  prediction 

of  concrete  dielectric  properties  in  the  high  frequency  range. 

This  study  shows  that  the  seasonal  temperature  variation 
usually  encountered  in  field  has  negligible  effect  on  the  real  part  of 

concrete  dielectric  permittivity  but  has  a  significant  impact  on  the 

attenuation  coefficient  of  concrete,  because  the  water  dielectric 
conductivity  drops  significantly  with  temperature.  A  temperature 
drop  from  20°C  to  5°C  causes  a  25%  reduction  in  the  attenuation 

coefficient  of  concrete.  Such  a  significant  drop  in  attenuation  has 

indeed  been  observed  in  bridge  decks  surveyed  using  radar  in 
summer  as  well  as  in  winter.  Besides,  the  attenuations  predicted 

with  the  new  model  are  closer  to  the  observed  field  data  than  the 
predictions  made  in  previous  studies  which  used  the  dc  conductivity 
instead  of  dielectric  conductivity. 
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The  theoretical  models  presented  in  this  study  facilitate  a 
rational  prediction  of  dielectric  properties  of  concrete  based  on  the 
volumetric  fraction  and  properties  of  its  constituents.  Future 
research  should  focus  on  the  verification  of  these  models  by  means 
of  either  radar  field  data  from  actual  bridge  decks  or  small  scale 
laboratory  experiments. 
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coulomb 

=  amp/sec 

volt 

=  kg  •  m2  •  amp*^  • 

sec"  3 

farad 

=  amp2  •  sec4  •  kg-1 

•  m* 

henry 

=  kg  •  m2  •  amp'2  • 

sec-  2 

ohm 

=  kg  •  m2  •  amp'2  • 

sec-  3 

S(Siemens) 

mho  =  ohm'^ 
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APPENDIX  III 

NUMERICAL  PROCEDURE  FOR  OBTAINING 
ROOTS  OF  COMPTTX  EQUATIONS 

(a)  Obtaining  Real  and  Complex  Roots 

The  complex  mixture  law  equation  (like  the  one  in  Eq.  32) 
involves  the  unknown  complex  quantity  e  which  has  a  real  part  £' 
and  an  imaginary  part  e"  which  can  be  obtained  using  a  complex  root 
solver.  A  complex  equation  can  also  be  written  as  two  real 
equations  as  follows: 

Real  part  of  the  expression  =  0 

Imaginary  part  of  the  expression  =  0 
Thus,  a  subroutine  which  finds  roots  (e'  and  e")  of  a  set  of  non-linear 
equations  can  also  be  employed.  Most  of  the  standard  subroutines 
require  good  initial  guesses  for  the  solution  to  converge.  Simple 
closed  form  equations  like  the  CRIM  model  (see  Eq.  27)  can  be  used 
for  obtaining  such  initial  guesses  for  the  roots.  To  be  sure  that  the 
root  has  converged,  it  is  advisable  to  use  more  than  one  root  solver 
for  the  same  equation. 

(b)  Computer  Programs  and  Subroutines  (FORTRAN-T?) 

The  following  computer  programs  (all  with  ’Double  Precision') 


have  been  used: 
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(i)  M3»F0R  and  M4*F0R  -  Water  conductivity,  CRIM  and 
Continuous  Grain  Size  Models  (the  two  programs  use 
different  complex  root  solvers) 

(ii)  M5*FOR  -  Water  Conductivity,  CRIM  and  Discrete  Grain 
Size  Models 

(iii)  PLl*FOR  -  Program  for  plotting  curves 

These  programs  are  listed  in  the  following  pages.  The  root  solving 
programs  employ  the  IMSL  subroutines  ZSCNT  and  ZSPOW  (see 
attached  description)  to  solve  for  roots  of  a  set  of  non-linear 
equations.  Further  documentation  on  IMSL  library  of  subroutines 
is  available  in  MIT's  Micro  VAX  Computer  Center  (Room  1-249). 
The  plotting  programs  use  the  'Penplot'  library  subroutines  which 
are  also  available  in  the  above  mentioned  computer  center. 
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A************************************************************ 

*  PROGRAM  NAME  - >  M3 . FOR - >  NO  SCREEN  INPUT 

*  SENDS  ADDITIONAL  OUTPUT  IN  'M2. DAT'  FOR  PLOTTING 

*  THIS  PROGRAM  INCORPORATES  THEORETICAL  MODELS  FOR 

*  (A)  CONDUCTIVITY  OF  SALINE  WATER 

*  (B)  DIELECTRIC  PROPERTIES  OF  SALINE  WATER 

*  (C)  MIXTURE  LAWS  FOR  CONCRETE 

*  INPUT  IS  READ  FROM  TERMINAL 

*  OUTPUT  IS  SENT  TO  FILES  CALLED  'Ml. DAT'  AND  'M2. DAT' 

*  WHICH  SHOULD  BE  OPENED  BEFORE  RUNNING  THIS  PROGRAM 
************************************************************* 

*  THIS  PROGRAM  USES  'IMSL'  SUBROUTINE  '2SCNT'  TO 

*  EVALUATE  COMPLEX  ROOTS 

*********************»*******'’************>^**’************»^*** 
A,B  -  FUNCTIONS  WHICH  DEPEND  ON  'SSW'  AND  'T' 

C  -  FUNCTION  OF  'ESW,EPSW,  F  AND  TSW' 

Cl  -  FUNCTION  OF  'DEL'  AND  'NSW' 

AA,BB,CC,DD  -  FUNCTIONS  OF  'ETA,  ESOL  AND  EAIR' 

ATT  -  ATTENUATION  PARAMETER 
DEL  -  25.0  -  T 

EO  -  DIELECTRIC  PERMITTIVITY  OF  FREE  SPACE 
EPSW  -  PARAMETER  AT  INFINITE  FREQUENCY 
FOR  BOTH  PURE  OR  SALINE  WATER 
EPW  =  PARAMETER  FOR  PURE  WATER  AT  ZERO  FREQUENCY 
EAIR  «  DIELECTRIC  PERMITTIVITY  OF  AIR 
EPS  =  COMPLEX  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
EPSI  -  IMAGINARY  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
EPSR  -  REAL  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
ESOL  -  DIELECTRIC  PERMITTIVITY  OF  SOLID 
ESW  =  PARAMETER  FOR  SALINE  WATER  AT  ZERO  FREQUNCY 
ETA  =  FUNCTION  OF  'P'  AND  'S' 

EW  «  COMPLEX  DIELECTRIC  PERMITTIVITY  FOR  WATER 

EWI  -  IMAGINARY  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  WATER 

EWI1,EWI2  -  TWO  COMPONENTS  OF  EWI 

EWR  -  REAL  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  WATER 
F  =  FREQUENCY  (GHz) 

ri,F2,F3,F4,F(l) ,F(2)  =  FUNCTIONS  DEFINING  COMPLEX 

EQUATION  FOR  THE  SPHERICAL  MODEL 
FCN  -  EXTERNAL  FUNCTION  FOR  IMSL  SUBROUTINE  'ZSCNT' 

FNORM  -  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCNT' 

I  =  AN  INTEGER  NUMBER 

lER  =  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCNT' 

IIM  =  IMAGINARY  NUMBER  =  (0. 0,1.0) 

ITMAX  =  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCNT' 

J  -  AN  INTEGER  NUMBER 
K  -  COMPLEX  WAVE  NUMBER 

Ml  -  FUNCTION  OF  ' AA , BB , CC , EAI R  AND  ESOL' 

MUO  -  PERMEABILITY  OF  FREE  SPACE 
N  *  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCNT' 

NSIG  -  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCN''-' 

NSW  -  NORMALITY  OF  SALINE  WATER 
P  *  POROSITY  OF  MIXTURE 

PAR  =  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCNT' 

PI  -  CONSTANT  =  3.141592654 
S  »  DEGREE  OF  SATURATION 


nono  nno  nn  no  on  o  n n n n o o o n n o n n 
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SI  -  FUNCTION  WHICH  DEPENDS  ON  'SSW'  AND  'T' 

SIG25  -  CONDUCTIVITY  AT  25  C 

SIGT  -  CONDUCTIVITY  AT  TEMPERATURE  'T' 

SSW  -  SALT  CONTENT  OF  WATER  [gm/kg] 

T  -  WATER  TEMPERATURE 

TPW  »  RELAXATION  TIME  FOR  PURE  WATER 
TSW  -  RELAXATION  TIME  FOR  SALINE  WATER 
W  -  ANGULAR  FREQUENCY  =  2.0*PI*F 
WK  »  PARAMETER  FOR  IMSL  SUBROUTINE  '2SCNT' 

X  »  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSCNT' 
********★****★*****★**********»***’'*******»'*’>******♦****♦  »*** 

234567890 

INTEGER  I , J,NSIG,N,ITMAX,IER 

REAL* 8  A,B,C,C1 , AA, BB , CC , DD , ATT , DEL 

REAL* 8  EAIR,E0 , EPSR , EPSI , EPSW , EPW 

REAL*8  ESOL , ESW , ETA , EWI , EWI 1 , EWI 2 , EWR 

REAL* 8  F,M1 ,MU0 ,NSW,P,PI ,S,SI,SIG25,SIGT 

REAL*8  SSW, T, TPW, TSW, W 

REAL*8  PAR{ 2 ) ,X( 2 ) , FNORM,WK( 42 ) 

COMPLEX*16  EW,EPS,IIM,K 
EXTERNAL  FCN 

COMMON/Ul/P , S , ESOL , EAIR , EW 
OPEN(UNIT-4 , FILE='M1 .DAT' , STATUS= ' OLD ' ) 

OPEN ( UNI T- 3, FI LE=' M2. DAT' , STATUS= ' OLD ' ) 

EO  -  8.854E-12 
PI  -  3.141592654 
MUO  =  PI*4.0E-7 

READ  INPUT  PARAMETERS 

★  ★★★Hr**************** 

S  *  0.0 
DO  20  1*1,6 
F  *  0.2 
DO  30  J=l,15 

PRINT  *, 'INPUT  VALUE  OF  SSW ( ppt ) , T { C )  AND  F(GH2)' 

READ (5,*)  SSW,T,F 
SSW-80. 0 
T*20 . 0 
F=F*1 . 0E9 

PRINT  *, 'INPUT  VALUE  OF  ?,S' 

READ( 5,* )P,S 
P=0.10 

INPUT  DATA 

EAIR-1 . 0 
ESOL-5 . 0 

COMPUTE  WATER  CONDUCTIVITY  (MODEL  FOR  'SSW'  LESS 
THAN  35.0  grams/liter) 

DEL  *  25.0  -  T 

SI  *  2.033E-2  +  1.266E-4*DEL  +  2 . 464E-6*DEL*DEL 

SI  -  SI  -  SSW* ( 1 . 849E-5  -  2.55lE-7*DEL  +  2 . 551E-8*DEL*DEL ) 


ononn  noon  o  o  nnooo 


-  71  - 


Page  3 


SI  -  DEL*SI 

SIG25  «  0.18252  -  1 . 4619E-3*SSW  +  2 . 093E-5*SSW*SSW 
S1G25  -  SSW*(S:G25  -  1 . 282E-7 * ( SSW* *  3 ) ) 

SIGT  -  SIG25*DEXP(-SI ) 


COMPUTE  REAL  AND  IMAGINARY  PARTS  OF  DIELECTRIC 
PEF.;;iTTIVITy  for  saline  water  (MODEL  FOR  'SSW' 


LESS  THAN  35.0  grams/liter) 


B  -  1.0  +  2 . 282E-5*T*SSW  -  7.638E-4*SSW 
B  «  B  -  7 . 760E-6*SSW*SSK  +  1 . 1 0 5E-8 * ( SSW* *  3 ) 
TPW  -  1.1109E-10  -  3.8/4E-12*T  +  6 . 938E-1 4 *T*T 
TPW  -  (TPW  -  5.096E-16*(T**3) )/(2.0*PI ) 

TSW  =  TPW*B 


A  -  1.0  +  1 . 613E-5*T*SSW  -  3.656E-3*SSW 
A  »  A  +  3 .21E-5*SSW*SSW  -  4 . 2 32E-7 * ( SSW* *  3 ) 
EPW  -  87.134  -  1.949E-1*T  -  1.276E-2*T*T 
EPW  -  EPW  +  2.491E-4*(T**3) 

ESW  =  EPW*A 


NSW  =  1.707E-2  +  1.205E-5*SSW  +  4 . 0 5 8E-9 * SSW* SSW 
NSW  =  NSW*SSW 
IF(SSW.LE.35.0)GO  TO  100 


COMPUTE  WATER  CONDUCTIVITY  (MODEL  FOR  'SSW' 
BETWEEN  35.0  AND  157.0  grams/liter) 

Cl  -  1.0  -  1.96E-2*DEL  +  8 . 0eE-5*DEL*DEL 
Cl  -  Cl  -  NSW*DEL* ( 3 . 02E-5  +  3.92E-5*DEL) 

Cl  =  Cl  -  NSW*NSW*DEL*( 1 .72E-5  -  6.58E-6*DEL) 
S1G25  =  10.39  -  2  37e*NSW  +  0 . 683*NSW*KSW 
SIG25  =  SIG25  -  0 . 135* (NSW**3 )  +  1 . OlE-2 * ( NSW* * 4 ) 
SIG25  =  NSW*SIG25 
SIGT  =  SIG25*C1 


COMPUTE  REAL  AND  IMAGINARY  PARTS  OF  DIELECTRIC 
PERMITTIVITY  FOR  SALINE  WATER  (MODEL  FOR  'SSW' 

BETWEEN  35.0  AND  157.0  crams/liter) 

B  «  1.0  +  0 . 146E-2*T*KSW  -  4.e96E-2*KSW 
B  =  B  -  2  .  97E-2*NSW*NSW  5 . 6  4 E- 3  *  (  NSW*  *  3  ) 

TSW  =  TPW*B 
C 

A  -  1.0  -  0.255*NSW  +  5 . 15E-2*KSW*NSW  -  6 . 6 9E- 3  * ( NSW* *  3 ) 
EPW  »  88.045  -  0.4147*T  +  6.295E-4*T*T  +  1  . 075E-5*  ( T*’'3  ) 
ESW  -  EPW*A 
C 

100  CONTINUE 

EPSW  -  4.9 
W  -  2.0*PI*F 

C  -  ( ESW- EPSW) /( 1 . 0  +  W*W-TSW*TSW) 

EWR  -  EPSW  +  C 
EWIl  «  W*TSW*C 
EWI2  =  SIGT/'W*E0) 
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EWI  -  EWIl  +  EWI2 
IIM  -  (0.0, 1.0) 

EW  -  EWR  +  IIM* EWI 


COMPLEX-REFRACTIVE  INDEX  METHOD  (CRIM) 


EPS  -  P*S*CDSQRT( EW)  +  P* ( 1 . 0- S ) * SQRT ( EAI R ) 
EPS  =  EPS  +  (1.0-P)*SQF.T(ESOL) 

EPS  -  EPS*EPS 
EPSR  »  DREAL(EPS) 

EPSI  -  DIMAG(EPS) 


K  -  W*CDSQRT( EPS*MUO*EO) 

ATT  -  DIMAG(K) 

OUTPUT  OF  NUMERICAL  RESULTS 


CONDUCTIVITY  MODEL 


F  =  F/1.0E9 
WRITE( 4,18) 

WRITE( 4,10) 

WRITE ( 4,11)  SSW,NSW,T, F, SIGT, EWR, EWIl , EWI 2 , EWI 
8  FORMAT(/4X, 'CONDUCTIVITY  MODEL') 

0  FORMAT ( 4X, ' SSW' , 5X, 'NSW' , 5X, ' T' ,5a, ' F ( GHz ) ' , 4X, 'SIGT' , 5X, ' EWR' , 

&4X, ' EWIl ' , 4X, ' EWI2 ' , 7X, ' EWI ' ) 

1  FORMAT{7F8.2,2E11.3) 

CRIM  MODEL 


WRITE( 4 , 19 ) 

WRITE( 4 , 12 ) 

WRITE (4, 13)  EAIR,ESOL,P,S,EPSR,EPSI ,ATT 
9  FORMAT(/4X, ' CRIM  MODEL') 

2  FORMAT ( 4X, ' EAIR ' , 4X, ' ESOL ' , 7X , ' P ' , 7a , ' S ' , 4X , ' EPSR' , 
&4X, 'EPSI' ,7X, 'ATT' ) 

3  FORMAT( 6F8 . 2 , Ell  .  3 ) 

SPHERICAL  GRAIN  SHAPE  MODEL 

INPUT  PARAMETERS  FOR  SUBROUTINE  'ZSCNT'  WHICH 
SOLVES  FOR  COMPLEX  ROOTS  OF  CON.PLEX  EQUATION 


N  »  2 
NSIG  -  3 
ITMAX  -  200 
PAR { 1 )  »  0.0 
PAR ( 2 )  -  0.0 

INITIAL  GUESS  FOR  ROOTS 


X(l)  -  EPSR 
X(2)  -  EPSI 
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SOLVE  FOR  COMPLEX  'EPS'  USING  SUBROUTINE  '2SCNT' 


CALL  2SCNT( FCN , NS IG , N , I TMAX , PAR , X , FNORM,WK, lER) 

EPS  =  X( 1 )  +  I1M*X( 2 ) 

K  »  W’'CDSQRT(  EPS*MU0*E0  ) 

ATT  »  DIMAG(K) 

OUTPUT  FOR  SPHERICAL  MODEL 


WRITE( 4,14) 

WRITE( 4 , 15 ) 

WRITE( 4 ,16)X(1),X(2) ,ATT 
WRITE( 4 , 17 ) 

4  FORMAT(/4X, 'SPHERICAL  MODEL') 

5  FORMAT (4X, 'EPSR' ,4X, 'EPSI' ,7X, 'ATT' ) 

6  FORMAT(2F8.2,E11.3) 

7  format!'*  ********************'//) 
WRITE(  3 , 2'?)F,X(T)  ,X(  2  )  ,ATT 

5  FORMAT ( 3 F8. 2, El 1.3) 

F  =  F  +  0.2 
0  CONTINUE 

S  -  S  +  0 . 2 
0  CONTINUE 


STOP 

END 

SUBROUTINE  — >  EQUATION  FOR  SPHERICAL  MODEL 


SUBROUTINE  FCN ( X , F , N , PAR ) 

COMPLEX* 16  EW, EPS ,Fl,F2,F3,r4,IIM 
REAL*8  AA,BB,CC,DD,ETA,M1 
REAL*8  P,S,ESOL,EAIR 
REAL*8  X( 2 )  , F( 2  )  ,PAR( 2 ) 

COMMON/Ul/P , S , ESOL , EAIR , EW 
C 

ETA  -  P*(1.0-S)/(1.0-P) 

AA  -  -2.0*(1.0-ETA) 

BB  =  2.0*ESOL  -  EAIR  +  ETA* ( 2 . 0*EAIR-ESOL) 
CC  -  ESOL*EAIR* ( 1 . 0+ETA) 

DD  *  BB*BB-4 . 0*AA*CC 
DD  -  SQRT(DD) 

Ml  -  2.0*ESOL  +  EAIR  +  ETA* ( 2 . 0 * EAI R+ESOL ) 
Ml  -  -M1/(2.0*DD) 

C 

IIM  -  (0.0, 1.0) 

EPS  -  X(l)  +  IIM*X(2) 

FI  -  AA*EPS*EPS  +  BB*EPS  +  CC 

F2  -  AA*EW*EW  +  BB*EW  +  CC 

FI  -  ( (EW/EPS )**(!. 0/3.0 ) )»CDSQRT(F1/F2) 
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C 

F3  -  2.0*AA*EPS  +  BB  -  DD 

F4  -  2.0*AA*EPS  +  BB  +  DD 

F3  -  F3  *  (2.0*AA*EW  +  BB  +  DD ) 

F4  »  F4  *  (2.0*AA*EW  +  BB  -  DD ) 

F3  =  (F3/F4)**M1 
C 

Fl  •=  F1*F3 

F(l)  -  DREAL(Fl)  -  S*P 
F(2)  -  DIMAG(Fl) 

C 

RETURN 

END 
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*  PROGRAM  NAME  - >  M4.F0R - >  NO  SCREEN  INPUT 

*  SENDS  ADDITIONAL  OUTPUT  IN  'M2. DAT'  FOR  PLOTTING 

*  THIS  PROGRAM  INCORPORATES  THEORETICAL  MODELS  FOR 

*  (A)  CONDUCTIVITY  OF  SALINE  WATER 

*  (B)  DIELECTRIC  PROPERTIES  OF  SALINE  WATER 

*  (C)  MIXTURE  LAWS  FOR  CONCRETE 

*  INPUT  IS  READ  FROM  TERMINAL 

*  OUTPUT  IS  SENT  TO  FILES  CALLED  'Ml. DAT'  AND  'M2. DAT' 

*  WHICH  SHOULD  BE  OPENED  BEFORE  RUNNING  THIS  PROGRAM 
************************************************************* 

*  THIS  PROGRAM  USES  'IMSL'  SUBROUTINE  'ZSPOW'  TO 

*  EVALUATE  COMPLEX  ROOTS 

************************************************************* 
A,B  -  FUNCTIONS  WHICH  DEPEND  ON  'SSW'  AND  'T' 

C  »  FUNCTION  OF  'ESW,EPSW,  F  AND  TSW' 

Cl  -  FUNCTION  OF  'DEL'  AND  'NSW' 

AA,BB,CC,DD  -  FUNCTIONS  OF  'ETA,  ESOL  AND  EAIR' 

ATT  -  ATTENUATION  PARAMETER 
DEL  -  25.0  -  T 

EO  -  DIELECTRIC  PERMITTIVITY  OF  FREE  SPACE 
EPSW  =  PARAMETER  AT  INFINITE  FREQUENCY 
FOR  BOTH  PURE  OR  SALINE  WATER 
EPW  =  PARAMETER  FOR  PURE  WATER  AT  ZERO  FREQUENCY 
EAIR  -  DIELECTRIC  PERMITTIVITY  OF  AIR 
EPS  -  COMPLEX  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
EPSI  =  IMAGINARY  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
EPSR  =  REAL  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
ESOL  =  DIELECTRIC  PERMITTIVITY  OF  SOLID 
ESW  =  PARAMETER  FOR  SALINE  WATER  AT  ZERO  FREQUNCY 
ETA  =  FUNCTION  OF  'P'  AND  'S' 

EW  -  COMPLEX  DIELECTRIC  PERMITTIVITY  FOR  WATER 

EWI  -  IMAGINARY  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  WATER 

EWI1,EWI2  «  TWO  COMPONENTS  OF  EWI 

EWR  -  REAL  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  WATER 
F  -  FREQUENCY  (GHz) 

Fl,F2,F3,F4,F(l) ,F(2)  =  FUNCTIONS  DEFINING  COMPLEX 

EQUATION  FOR  THE  SPHERICAL  MODEL 
FCN  -  EXTERNAL  FUNCTION  FOR  IMSL  SUBROUTINE  'ZSPOW' 

FNORM  -  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

I  -  AN  INTEGER  NUMBER 

lER  *  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

IIM  =  IMAGINARY  NUMBER  =  (0,0, 1.0) 

ITMAX  =  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

J  *  AN  INTEGER  NUMBER 
K  -  COMPLEX  WAVE  NUMBER 

Ml  »  FUNCTION  OF  '  AA ,  BB  ,  CC  ,  E^.  I R  AND  ESOL' 

MUO  -  PERMEABILITY  OF  FREE  SPACE 
N  »  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

NSIG  -  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

NSW  -  NORMALITY  OF  SALINE  WATER 
P  -  POROSITY  OF  MIXTURE 

PAR  -  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

PI  -  CONSTANT  -  3.141592654 
S  -  DEGREE  OF  SATURATION 
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SI  -  FUNCTION  WHICH  DEPENDS  ON  'SSW'  AND  'T' 

SIG25  -  CONDUCTIVITY  AT  25  C 

SIGT  «  CONDUCTIVITY  AT  TEMPERATURE  'T' 

SSW  -  SALT  CONTENT  OF  WATER  [gm/kg] 

T  -  WATER  TEMPERATURE 

TPW  -  RELAXATION  TIME  FOR  PURE  WATER 
TSW  *  RELAXATION  TIME  FOR  SALINE  WATER 
W  -  ANGULAR  FREQUENCY  -  2.0*PI*F 
WK  »  PARAMETER  FOR  IMSL  SUBROUTINE  '2SP0W' 

X  -  PARAMETER  FOR  IMSL  SUBROUTINE  'ZSPOW' 

★  ***l>************»*lk***************************************** 

567890 

INTEGER  I,J,NSIG,N,ITKAX,IER 

REAL* 8  A,B,C,C1 , AA, BB , CC , DD , ATT , DEL 

REAL* 8  EAIR,E0 ,EPSR,EPSI ,EPSW,EPW 

REAL* 8  ESOL , ESW , ETA , EWI , EWI 1 , EWI 2 , EWR 

REAL* 8  F,M1,MJ0,NSW,P,PI,S,SI, SIG25,EIGT 

REAL*8  SSW, T, TPW, TSW, W 

REAL*8  PAR( 2 ) ,X( 2 ) ,FNORM,WK( 21 ) 

COMPLEX*16  EW,EPS,IIM,K 
EXTERNAL  FCN 

COMMON/Ul /P , S , ESOL , EAIR , EW 

OPEN (UNIT- 4, FILE- 'Ml. DAT' , STATUS- ' OLD ' ) 

OPEN(UNIT=3,FILE='M2.DAT' , STATUS- ' OLD ' ) 

EO  -  8.854E-12 
PI  =  3.141592654 
MUO  =  P1*4.0E-7 

READ  INPUT  PARAMETERS 

S  -  0.0 
DO  20  1-1,6 
F  -  0.2 
DO  30  J-1,15 

PRINT  *, 'INPUT  VALUE  OF  SSW ( ppt ) , T ( C )  AND  F{GH2)' 

READ(5,*)  SSW,T,F 
SSW-12. 0 
T-20.0 
F=F*1 . 0E9 

PRINT  *, 'INPUT  VALUE  OF  ?,S' 

READ( 5 , * ) P, S 
P-0.10 

INPUT  DATA 

EAIR-1 . 0 
ESOL-5.0 

COMPUTE  WATER  CONDUCTIVITY  (MODEL  FOR  'SSW'  LESS 
THAN  35.0  grams/liter ) 

DEL  -  25.0  -  T 

SI  »  2.033E-2  +  1.266E-4*DEL  +  2 . 464E-6*DEL*DEL 

SI  -  SI  -  SSW* ( 1 . 849E-5  -  2.55lE-7*DEL  +  2 . 551E-8*DEL*DEL ) 
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SI  -  DEL*SI 

SIG25  -  0.18252  -  1 . 4619E-3*SSW  +  2 . 093E-5*SSW*SSW 
SIG25  -  SSW*(SIG25  -  1 . 282E-7 * ( SSW* * 3 ) ) 

SIGT  -  SIG25*DEXP(-SI ) 

C 

C  COMPUTE  REAL  AND  IMAGINARY  PARTS  OF  DIELECTRIC 

C  PERMITTIVITY  FOR  SALINE  WATER  (MODEL  FOR  'SSW' 

C  LESS  THAN  35.0  grams/liter) 

Q  lfc********^******'*******«*K***'^r************'#t** 

B  -  1.0  +  2.282E-5*T*SSW  -  7.638E-4*SSW 
B  -  B  -  7 .760E-6*SSW*SSW  +  1 . 105E-8* ( SSW**3 ) 

TPW  -  1.1109E-10  -  3.824E-12*T  +  6 . 938E-14^T*T 
TPW  -  (TPW  -  5.096E-16*(T**3) )/(2.0*PI) 

TSW  »  TPW*B 
C 

A  «  1.0  +  1 . 613E-5*T*SSW  -  3.656E-3*SSW 
A  -  A  +  3 . 21E-5*SSW*SSW  -  4 . 2 32E-7 * ( SSW* *  3 ) 

EPW  -  87.134  -  1.949E-1*T  -  1.276E-2*T*T 
EPW  -  EPW  +  2 . 491E-4* ( T**3 ) 

ESW  -  EPW*A 
C 

NSW  -  1.707E-2  +  1.205E-5*SSW  +  4 . 058E-9*SSW*SSW 
NSW  -  NSW*SSW 
IF( SSW.LE. 35 . 0 )GO  TO  100 
C 

C  COMPUTE  WATER  CONDUCTIVITY  (MODEL  FOR  'SSW' 

C  BETWEEN  35.0  AND  157.0  grams/liter) 

Cl  -  1.0  -  1.96E-2*DEL  +  8 . 08E-5*DEL*DEL 
Cl  -  Cl  -  NSW*DEL*(  3.02E-5  +  3 . 92E~5’'DEL  ) 

Cl  -  Cl  -  NSW*NSW*DEL* ( 1 .72E-5  -  6.58E-6*DEL) 

SlG2b  -  10.39  -  2.378*NSW  +  0 . 683*NSW*NSW 
S1G25  -  SIG25  -  0 . 1 35* ( NSW** 3 )  +  1 . 0 lE-2 * ( NSW* * 4 ) 

S1G25  «  NSW*SIG25 
SIGT  =  SIG25*C1 
C 

C  COMPUTE  REAL  AND  IMAGINARY  PARTS  OF  DIELECTRIC 

C  PERMITTIVITY  FOR  SALINE  WATER  (MODEL  FOR  'SSW' 

C  BETWEEN  35.0  AND  157.0  grams/l i te r ) 

B  =  1.0  +  0 . 146E-2*T*NSW  -  4.696E-2*NSW 
B  =  B  -  2 . 97E-2*NSW*NSW  +  5 . 6 4 E- 3  * ( NSW* *  3 ) 

TSW  =  TPW*B 
C 

A  -  1.0  -  0.255*NSW  +  5 .15E-2*NSW*NSW  -  6 . 89E-3* ( NSW**3 ) 
EPW  -  88.045  -  0.4147*T  +  6.295E-4*T*T  +  1 . 075E-5* ( T**3 ) 
ESW  -  EPW*A 

r 

w 

100  CONTINUE 

EPSW  ■  4.9 
W  «  2.0*PI*F 

C  ■  ( ESW-EPSW)/( 1 . 0  +  W*W*TSW*TSW) 

EWR  -  EPSW  +  C 
EWIl  -  W*TSW*C 
EWI2  »  SIGT/(W*E0) 
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EWI  =  EWIl  +  EWI2 
IIM  -  (0.0, 1.0) 

EW  -  EWR  +  IIM*EWI 


COMPLEX-REFRACTIVE  INDEX  METHOD  (CRIK) 

EPS  -  P*S*CDSQRT( EW)  +  P* ( 1 . 0-S ' *SQRT ( EAIR ) 
EPS  =  EPS  +  { 1 . 0-P ) *SQRT( ESOL) 

EPS  »  EPS*EPS 
EPSR  -  DREAL(EPS) 

EPSI  -  DIMAG(EPS) 

K  -  W*CDSQRT( EPS*MU0^E0 ) 

ATT  »  DIMAG(K) 


OUTPUT  OF  NUMERICAL  RESULTS 


CONDUCTIVITY  MODEL 


F  -  F/1.0E9 
WRITE( 4 ,18 ) 

WRITE ( 4,10) 

WRITE( 4 , 11 )  SSW,NSW, T, F, SIGT, EWR, EWI 1 , EWI2 , EWI 
8  FORMAT(/4X, ' CONDUCTIVITY  MODEL') 

0  FORMAT ( 4X, ' SSW' , 5X , 'NSW' , 5X , ' T ' , 5X , ’ F ( GKz ) ' , 4X , ' SI GT ' , 5X , 

&4X, ' EWIl ' , 4X, ' EWI2 ' , 7X, ' EWI ' ) 

1  FORMAT! 7F8 . 2 , 2E11 , 3 ) 

CRIM  MODEL 


WRITE! 4, 19) 

WRITE! 4 , 12 ) 

WRITE! 4,13)  EAIR,ESOL,P,S,EPSR,EPSI ,ATT 
9  FORMAT(/4X, 'CRIM  MODEL') 

2  FORMAT! 4X, 'EAIR' ,4X, 'ESOL' , 7X , ' P ' , 7X , ' S ' , 4X , 'EPSR' , 
&4X, 'EPSI ' ,7X, 'ATT' ) 

3  FORMAT! 6F0 , 2 , Ell . 3 ) 

SPHERICAL  GRAIN  SHAPE  MODEL 


INPUT  PARAMETERS  FOR  SUBROUTINE  'ISPOW'  WHICH 
SOLVES  FOR  COMPLEX  ROOTS  OF  COMPLEX  EQUATION 


N  =  2 
NSTG  =  3 
ITMAX  -  200 
PAR ! 1 )  «  0.0 
PAR! 2)  »  0.0 

INITIAL  GUESS  FOR  ROOTS 


X!l)  *  EPSR 
X!2)  =  EPSI 


'  EWR '  , 
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SOLVE  FOR  COMPLEX  'EPS'  USING  SUBROUTINE  '2SP0W' 


CALL  ZSPOW{ FCN,NSIG,N, ITMAX,PAR,X, FNORM,WK, lER) 

EPS  -  X(l)  +  IIM*X(2) 

K  -  W*CDSQRT( EPS*MU0*E0 ) 

ATT  -  DIMAG(K) 

OUTPUT  FOR  SPHERICAL  MODEL 


WRITE( 4 , 14 ) 

WRITE( 4 , 15 ) 

WRITE ( 4 ,16)X(1),X(2) ,ATT 
WRITE( 4 ,17 ) 

4  FORMAT (/4X, 'SPHERICAL  MODEL' ) 

5  FORMAT (4X, 'EPSR' ,4X, 'EPSI' ,7X, 'ATT' ) 

6  FORMAT(2F8.2,E11.3) 

7  FORMAT{'* 

WRITE( 3,2^)F,X( 1 ) ,X( 2 ) ,ATT 
5  FORMAT(3F8.2,E11.3) 

F  =  F  +  0.2 
0  CONTINUE 

S  -  S  +  0.2 
0  CONTINUE 


STOP 

END 

SUBROUTINE  — >  EQUATION  FOR  SPHERICAL  MODEL 
******************************************** 


SUBROUTINE  FCN ( X , F , N , PAR ) 

COMPLEX* 16  EW,EPS ,F1,F2,F3,F4,IIM 
REAL*8  AA,BB,CC,DD,ETA,M1 
REAL*8  P,S,ESOL,EAIR 
REAL*8  X( 2 ) , F( 2 ) , PAR{ 2 ) 

COMMON/Ul/P , S , ESOL , EAIR , EW 
C 

ETA  »  P*(1.0-S)/(1.0-P) 

AA  =  -2.0*(1.0+ETA) 

BB  -=  2.0*ESOL  -  EAIR  +  ETA*  ( 2 . 0*EAIR-ESOL) 
CC  -  ESOL*EAIR* ( 1 . 0+ETA) 

DD  -  BB*BB-4 . 0*AA*CC 
DD  -  SQRT(DD) 

Ml  -  2.0*ESOL  +  EAIR  +  ETA* ( 2 . 0 * EAI R+ESOL ) 
Ml  -  -M1/(2.0*DD) 

C 

IIM  -  (0.0, 1.0) 

EPS  «  X(l)  +  IIM*X(2) 

FI  »  AA*EPS*EPS  +  BB*EPS  +  CC 

F2  ■  AA*EW*EW  +  BB*EW  +  CC 

FI  -  ( (EW/EPS)**(1.0/3.0) )  CDSQRT(F1/F2) 
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C 

F3  -  2.0*AA*EPS  +  BB  -  DD 

F4  -  2.0*AA*EPS  +  BB  +  DD 

F3  -  F3  *  (2.0*AA*EW  +  BB  +  DD ) 

F4  -  F4  *  (2.0*AA*EW  +  BB  -  DD ) 

F3  -  (F3/F4)**M1 
C 

FI  -  Fl*F3 

F(l)  •=  DREAL(Fl)  -  S*P 
F(2)  -  DIMAG(Fl) 

C 

RETURN 

END 
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*  PROGRAM  NAME  - >  M5.F0R - >  DISCRETE  GRAIN  SIZES  MODEL 

*  SENDS  ADDITIONAL  OUTPUT  IN  'M2. DAT'  FOR  PLOTTING 

*  AND  IN  'M3. DAT'  FOR  CHECKING  THE  NUMERICAL  SOLUTION  STEPS 

*  THIS  PROGRAM  INCORPORATES  THEORETICAL  MODELS  FOR 

*  (A)  CONDUCTIVITY  OF  SALINE  WATER 

*  (B)  DIELECTRIC  PROPERTIES  OF  SALINE  WATER 

*  (C)  MIXTURE  LAWS  FOR  CONCRETE 

*  INPUT  IS  READ  FROM  TERMINAL 

*  OUTPUT  IS  SENT  TO  FILES  CALLED  ' Ml . DAT ' , ' M2 . DAT '  AND 

*  'M3. DAT'  WHICH  SHOULD  BE  OPENED  BEFORE  RUNNING  THIS  PROGRAM 
******************★****<>********************************★**** 

A,B  -  FUNCTIONS  WHICH  DEPEND  ON  'SSW'  AND  'T' 

AN(I)  -  FINAL  PROPORTION  OF  EACH  SIZE  COMPONENT 
C  -  FUNCTION  OF  'ESW,EPSW,  F  AND  TSW' 

Cl  -  FUNCTION  OF  'DEL'  AND  'NSW' 

AA,BB,CC  «  FUNCTIONS  OF  'PN(I)'  AND  'SIGC(l)' 

ATT  -  ATTENUATION  PARAMETER 
DEL  -  25.0  -  T 

EO  -  DIELECTRIC  PERMITTIVITY  OF  FREE  SPACE 
EPSW  -  PARAMETER  AT  INFINITE  FREQUENCY 
FOR  BOTH  PURE  OR  SALINE  WATER 
EPW  -  PARAMETER  FOR  PURE  WATER  AT  ZERO  FREQUENCY 
EAIR  =  DIELECTRIC  PERMITTIVITY  OF  AIR 
EPS  *  COMPLEX  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 
EPSCRM  -  SAME  AS  'EPSR'  FOR  S  <  0,01 

EPSI  «  IMAGINARY  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 

EPSR  =  REAL  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  MIXTURE 

ESOL  =  DIELECTRIC  PERMITTIVITY  OF  SOLID 

ESW  =  PARAMETER  FOR  SALINE  WATER  AT  ZERO  FREQUNCY 

EW  -  COMPLEX  DIELECTRIC  PERMITTIVITY  FOR  WATER 

EWI  «  IMAGINARY  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  WATER 

EWI1,EWI2  -  TWO  COMPONENTS  OF  EWI 

EWR  »=  REAL  PART  OF  DIELECTRIC  PERMITTIVITY  FOR  WATER 
F  =  FREQUENCY  (GHz) 

I,  I I  -  AN  INTEGER  NUMBER 

IIM  -  IMAGINARY  NUMBER  =  (0.0, 1.0) 

J, JJ,KK  -  AN  INTEGER  NUMBER 
K  «  COMPLEX  WAVE  NUMBER 

MUO  =  PERMEABILITY  OF  FREE  SPACE 
NSW  »  NORMALITY  OF  SALINE  WATER 
P  »  POROSITY  OF  MIXTURE 

PN(I)  =  PROPORTION  OF  INCLUSION  AT  EACH  STEP 
PI  =  CONSTANT  =  3.1415S2654 
S  «  DEGREE  OF  SATURATION 

SI  »  FUNCTION  WHICH  DEPENDS  ON  'SSW'  AND  'T' 

SIG25  -  CONDUCTIVITY  OF  WATER  AT  25  C 

SIGT  -  CONDUCTIVITY  OF  WATER  AT  TEMPERATURE  'T' 

SIGC(I)  -  CONDUCTIVITY  OF  INCLUSION  AT  EACH  STEP 
SIGTC  -  COMPLEX  CONDUCTIVITY  OF  THE  MIXTURE 
SIGRC  -  REAL  PART  OF  COMPLEX  CONDUCTIVITY 
SIGIC  =  IMAGINARY  PART  OF  COMPLEX  CONDUCTIVITY 
SSW  *  SALT  CONTENT  OF  WATER  (gm/)<g] 

T  »  WATER  TEMPERATURE 

TPW  -  RELAXATION  TIME  FOR  PURE  WATER 


nnoon  nonn  onono  no  no  o  noon 


-  82  - 


Page  2 


TSW  =  RELAXATION  TIME  FOR  SALINE  WATER 
W  -  ANGULAR  FREQUENCY  -  2.0*PI*F 

234567890 

INTEGER  I,J,II,JJ,KK 
REAL*8  A,B,C,C1,ATT,AN(11),DEL 
REAL*8  EAIR,EO,EPSCRM,EPSR,EPSI ,EPSW,EPW 
REAL  *  8  ESOL , ESW , EWI , EWI 1 , EWI 2 , EWR 
REAL*8  F,MU0,NSW,P,PI ,PN( 11 ) ,S,SI,SIG25 
REAL*8  SIGIC,SIGRC,SIGT,SSW,T,TPW,TSW,W 
COMPLEX*16  AA,BB,CC,EW,EPS,IIM,K 
COMPLEX*16  SIGC{ 11 ) ,SIGTC 

OPEN(UNIT-4 , FILE*'M1 .DAT' , STATUS= ' OLD ' ) 

OPEN ( UNI T-3, FILE- 'M2. DAT' , STATUS- ' OLD ' ) 
OPEN{UNIT-2,FILE-'M3.DAT' , STATUS- ' OLD ' ) 

EO  -  8.854E-12 
PI  -  3.141592654 
MUO  -  PI*4.0E-7 

READ  INPUT  PARAMETERS 

iritifitif-kifitifitirififititifitififific 

S  —  0.0 
DO  20  1=1,6 
F  -  0.2 
DO  30  J-1,15 

PRINT  *, 'INPUT  VALUE  OF  SSW ( ppt ) , T ( C )  AND  F(GHZ)' 
READ(5,*)  SSW,T,F 
SSW-52.0 
T-20 . 0 
F-F*1.0E9 

PRINT  *, 'INPUT  VALUE  OF  P , S ' 

READ( 5, * )P,S 
P-0.15 

INPUT  DATA 

EAIR-1 . 0 
ESOL-5.0 

COMPUTE  WATER  CONDUCTIVITY  (MODEL  FOR  'SSW'  LESS 
THAN  35.0  grams/liter ) 

DEL  -  25.0  -  T 

SI  »  2.033E-2  +  1.266E-4*DEL  +  2 . 464E-6*DEL*DEL 

SI  -  SI  -  SSW* ( 1 .849E-5  -  2.551E-7*DEL  +  2 . 551E-8*DEL*DEL ) 

SI  -  DEL*SI 

SIG25  -  0.18252  -  1 . 4619E-3'^SSW  +  2 . 093E-5*SSW*SSW 
SIG25  -  SSW*(SIG25  -  1 . 2B2E-7 * ( SSW* * 3 ) ) 

SIGT  -  S1G25*DEXP(-SI ) 

COMPUTE  REAL  AND  IMAGINARY  PARTS  OF  DIELECTRIC 
PERMITTIVITY  FOR  SALINE  WATER  (MODEL  FOR  'SSW' 

LESS  THAN  35.0  grams/liter) 
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B  -  1.0  +  2 . 282E-5*T*SSW  -  7.638E-4*SSW 
B  -  B  -  7.760E-6*SSW*SSW  +  1 . 1 0 5E-8* ( SSW* * 3 ) 

TPW  -  1.1109E-10  -  3.824E-12*T  +  6 . 9 38E-1 4 *T*T 
TPW  -  (TPW  -  5.096E-16*(T**3) )/(2.0*PI ) 

TSW  =  TPW*B 
C 

A  -  1.0  +  1 . 613E-5*T*SSW  -  3.656E-3*SSW 
A  -  A  +  3.21E-5*SSW*SSW  -  4 , 2 32E-7 * ( SSW* * 3 ) 

EPW  =  87.134  -  1.949E-1*T  -  1.276E-2*T*T 
EPW  =  EPW  +  2.491E-4*(T**3) 

ESW  «  EPW*A 
C 

NSW  -  1.707E-2  +  1.205E-5*SSW  +  4 . 058E-9*SSW*SSW 
NSW  «=  NSW*SSW 
IF(SSW.LE.35.0)GO  to  100 

COMPUTE  WATER  CONDUCTIVITY  (MODEL  FOR  'SSW' 

BETWEEN  35.0  AND  157,0  g r air.s/1  i te r  ) 

Cl  -  1.0  -  1.96E-2*DEL  +  8 . 08E-5*DEL*DEL 
Cl  -  Cl  -  NSW*DEL* ( 3 . 02E-5  +  3.92E-5*DEL) 

Cl  -  Cl  -  NSW*NSW*DEL*(1.72E-5  -  6.58E-6*DEL) 

SIG25  -  10.39  -  2.378*NSW  +  0 . 6 8 3 *NSW*NSW 
SIG25  =  SIG25  -  C . 1 3 5* ( NSW* *  3 )  +  1 . OlE-2* ( NSW**4 ) 

SIG25  =  NSW*SIG25 
SIGT  =  SIG25*C1 

COMPUTE  REAL  AND  IMAGINARY  PARTS  OF  DIELECTRIC 
PERMITTIVITY  FOR  SALINE  WATER  (MODEL  FOR  'SSW' 

BETWEEN  35.0  AND  157.0  gr ams/l i te r ) 

B  -  1.0  +  0 . 146E-2*T*NSW  -  4.896E-2*NSW 
B  -  B  -  2 . 97E-2*NSW*NSW  +  5 . 64E-3* ( NSW* *3 ) 

TSW  =  TPW*B 

A  =  1.0  -  0.255*NSW  +  5 . 15E-2*NSW*NSW  -  6 . 89E- 3 * ( NSW** 3 ) 
EPW  =  88.045  -  0.4147*T  +  6.295E-4*T*T  +  1 . 075E-5* ( T**3 ) 
ESW  =  EPW*A 

00  CONTINUE 

EPSW  =  4.9 
W  =  2.0*PI*F 

C  =  ( ESW-EPSW)/( 1 . 0  +  W*W*TSW*TSK) 

EWR  »  EPSW  +  C 
EWIl  -  W*TSW*C 
EWI2  -  SIGT/(W*E0) 

EVTT  -  EWIl  +  EWI2 
IIM  -  (0.0, 1.0) 

EW  -  EWR  +  IIM*EWI 

COMPLEX-REFRACTIVE  INDEX  METHOD  (CRIM) 

EPS  «  P*S*CDSQRT(  EW)  +  P*  (  1 . 0-S  )  *S(2RT (  FAIR  ) 

EPS  *  EPS  +  ( 1 . 0-P) *SQRT( ESOL) 

EPS  *  EPS*EPS 
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EPSR  •  DREAL(EPS) 

EPSI  -  DIMAG(EPS) 

IF(  S.LT.  0 . 01  )EPSCRM>=EPSR 


K  -  W*CDSQRT( EPS*MU0*E0 ) 

ATT  =  DIMAG ( K ) 

OUTPUT  OF  NUMERICAL  EESULCS 


CONDUCTIVITY  MODEL 


F  -  F/1.0E9 
WRITE! 4 ,18 ) 

WRITE! 4,10) 

WRITE! 4,11)  SSW,NSW,T,F, SIGT, EWR, EWIl , EWI2 , EWI 
FORMAT!/4X, ' CONDUCTIVITY  MODEL') 

FORMAT! 4X, ' SSW'  , 5X, 'NSW'  , 5X , ' T '  , 5X ,  ' F  !  GHz )  '  , 4X ,  ' S I GT '  , 5X , ' EWR 
&4X, 'EWIl' ,4X, 'EWI 2' ,7X, 'EWI' ) 

FORMAT! 7F8 . 2 , 2E11 . 3 ) 

CRIM  MODEL 


WRITE! 4,19) 

WRITE!4,12) 

WRITE! 4, 13)  EAIR,ESOL,P,S, EPSR, EPSI ,ATT 
FORMAT!/4X, 'CRIM  MODEL') 

FORMAT! 4X, ' EAIR' , 4X, ' ESOL' ,7X,'P',7X,'S',4X,'EPSR', 
&4X, 'EPSI ' ,7X, 'ATT' ) 

FORMAT !6F8. 2, El 1.3) 

DISCRETE  GRAIN  SIZE  MODEL 


INPUT  VALUES  FOR  FINAL  MIX  PROPORTIONS 


AN! 1 )  -  0 . 5/3 . 0 
AN! 2)  =  0. 5/3.0 
AN! 3)  =  0. 5/3.0 
AN!4)  =  ! 1 . 0-S ) *P/3 . 0 
AN!5)  =  !0.5  -  P)/2.0 
AN!6)  =  3.0*!0.5  -  P)/10.0 
AN!7)  =  2.0*!0.5  -  P)/10.0 
AN! 8)  =  ! 1 . 0-S ) *P/3 . 0 
AN!9)  =  2.0*!1.0-S)*P/9.0 
AN!10)  -  2.0*!1.0-S)*P/18.0 
AN ! 1 1  )  =  P  *  S 

INPUT  VALUES  FOR  COMPLEX  DIELECTRIC 
PERMITTIVITY  FOR  EACH  COMPONENT 


SIGC!1)  -  5.0 
SIGC!2)  »  5.0 
SIGC!3)  -  5.0 
SIGC! 4  )  =  1 . 0 


SIGC(5)  -  5.0 
SIGC(6)  -  5.0 
SIGC(7)  «  5.0 
SIGC(8)  -  1.0 
SIGC(9)  *  1.0 
SIGC(IO)  -  1.0 
SIGC(ll)  -  EW 
C 

C  CONVERT  THE  ABOVE  VALUES  TO 

C  COMPLEX  CONDUCTIVITIES 

C  - 

DO  40  II»1,11 

SIGC(II)  -  -IIM*W*E0’'S1GC{  II  ) 

40  CONTINUE 

C 

C  COMPUTE  PN(I) 

C  - 

PN(1)  -  AN(1) 

PN(2)  -  AN(2)/(1.0  PN(1)) 

DO  50  II  «  3,11 
PN ( I I )  -  AN ( I I ) 

DO  50  JJ  «  1,11-1 

PN(II)  »=  PN(II)/(1.0  -  PN(JJ)) 

50  CONTINUE 

C 

C  OUTPUT  FOR  CHECKING  THE  RESULTS 

C  - 

WRITE(2,21) 

WRITE(2,22) (AN(II ) ,11=1,11) 

WP  T  TP  O  7  ^ 

WRITE(2!22)(PN(II) ,11-1,11) 

21  FORMAT(/4X, ' (AN( I ) ,1-1,11) ' ) 

22  FORMAT(6F10.4 ) 

23  FORMAT(/4X, ' (PN(I) ,1  =  1,11)  '  ) 

C 

C  COMPUTE  COMPLEX  CONDUCTIVITY 

C  OF  THE  MIXTURE  AT  EACH  STEP 

C  - 

SIGTC  -  SIGC(ll) 

DO  60  11-1,10 
JJ  -  11  -  II 
AA  =  2 . 0 

CC  -  -SIGTC*SIGC( JJ) 

BB  -  3.0*(SIGTC-SIGC( JJ) )*PN( JJ) 

BB  -  BB  +  SIGC(JJ)  -  2.0*SIGTC 
C 

prnTC  m  CDFOPTt  ~  4.'^*»i»*CC) 

SIGTC  -  -(BB  -  SIGTC)/(2.0*AA) 

IF(AN(11) .LT. 0.001. AND. JJ.EQ.lOjSIGTC  «  SIGC(IO) 
SIGRC  -  DREAL( SIGTC) 

SIGIC  -  DIMAG(SIGTC) 

C 

C  OUTPUT  FOR  CHECKING  THE  ALGORITHM 

C  - - - 

WRITE(2,24) 
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WRITE(2,26)II,SIGRC,SIGIC 

24  F0RMAT(/4X, 'II' ,6X, 'SIGRC' ,7X, 'SIGIC' ) 

26  F0RMAT(I5,2E12.3) 

C 

60  CONTINUE 

C 

EPS  -  -SIGTC/( IIM*W*E0 ) 

EPSR  -  DREAL(EPS) 

EPSI  -  DIMAG(EPS) 

C 

K  «  W*CDSQRT( EPS*MUO*EO ) 

ATT  -  D1MAG{K) 

C 

C  OUTPUT  FOR  DISCRETE  GRAIN  SIZE  MODEL 

C  - 

WRITE(4,14) 

WRITE(4,15) 

IF(S.LT.O.Ol) EPSR-EPSCRM 
IF(S.LT,0.01)EPSI-0.0 
IF(S.LT.0.01)ATT-0.0 
WRITE (4, 16) EPSR, EPS I, ATT 
WRITE( 4,17) 

14  FORMAT(/4X, 'DISCRETE  GRAIN  SIZE  MODEL') 

15  FORMAT (4X, 'EPSR' ,4X, 'EPSI' ,7X, 'ATT' ) 

16  F0RMAT(2F8.2,E11.3) 

17  F0RMAT('* 

WRITE(  3,2'5)F,EPSR,EPSI  ,ATT 

25  F0RMAT(3F8.2,E11.3) 

C 

F  -  F  +  0.2 
30  CONTINUE 

S  *  S  +  0.2 
20  CONTINUE 

C 

STOP 

END 
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C  PROGRAM  NAME — >  PLl , FOR 

C  THIS  PROGRAM  PLOTS  RESULTS  OF  DIELECTRIC  DISPERSION  CURVES  FROM 
C  FILES  M2. DAT.  IT  USES  THE  PENPLOT  SUBROUTINE  'QPICTR' 

C  WHICH  PASSES  A  SMOOTH  CURVE  THROUGH  GIVEN  POINTS.  THIS  PARTICULAR 
C  VERSION  PLOTS  DIELECTRIC  PROPERTIES  VS.  FREQUENCY. 

C  BOTH  X&Y  AXES  ARE  IN  NORMAL  SCALES. 

c 

C23456789 

REAL  F,EPSR,EPSI ,ATT 
REAL  XLO,  YLO,XHI ,  YHI  ,DUMr:Y 
REAL  XSCL(4) ,BALK(2,15) 

INTEGER  I, J,K,IND(2) ,ISCL 
CHARACTER* 50  XLAB.YLAB 
C 

OPEN  (UNIT- 3, FILE- 'M2. DAT' , STATUS- ' OLD ' ) 

C 

QiH^**iHk****iimi***ii*ilii*****ii*i>*i>**iHiitiHHt*****iHt*iHt*****************ii* 

C  READ  PARAMETER  TO  DETERMINE  WHICH  PLOT  IS  REQUIRED 
PRINT  *, 'INPUT  PLOT  PARAMETER  K  &  YLO  &  YHI ' 

READ(5,*)  K, YLO, YHI 

ENTER  PARAMETERS  FOR  MINIMUM  AND  MAXIMUM  LIMITS  OF  THE  CURVE 
XLO-0 . 0 
XHI-3. 0 
YLO-0. 0 
YHI-10. 0 

ENTER  PARAMETERS  FOR  SMOOTH  CURVE 
ISCL— 2 
XSCL( 1 )-XLO 
XSCL( 2)-XHl 
XSCL( 3 )-YLO 
XSCL( 4 )»YHI 

LABELS 

XLAB- ' Frequency  (GHz)' 

IF ( K . EQ. 1 ) YLAB* ' Permittivity  of  Concrete{  )' 

IF ( K . EQ . 2 ) YLAB- ' Loss  Factor  of  Concrete(  )' 

IF(  K.EQ. 3 ) YLAB-'Attenuation  Coeff.  of  Concrete(  )' 


READ  DATA  FROM  INPUT  FILES 

READ  DATA  FROM  M2. DAT 

FIRST  READ  HEADER 
READ(3,30)DUMMY 

DO  50  1-1,6 
DO  40  J-1,15 

READ ( 3 , 30 ,END-45)F,EPSR,EPSI ,ATT 
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30  F0RMAT(3F8.2,E11.3) 

ATT-ABS ( ATT ) 

BALK( 1 , J )-F 

I F ( K . EQ . 1 ) BALK (2,3) -EPSR 
IF(K.EQ.2)BALK(2, J)-EPSI 
I F ( K . EQ . 3 ) BALK ( 2 , J ) -ATT 
40  CONTINUE 

45  1ND(1)-15 

1ND(2)-IND(1) 

PLOT  SMOOTH  CURVES  FOR  THE  GIVEN  DATA 


CALL  0PICTR(BALK,2,IND(1) ,Qy(2) ,QX(1) , QXSCL ( XSCL ) , 
&QISCL(1SCL), QMOVE (Cl), QLABEL (14), QXLAB ( XLAB ) , QYLAB ( YLAB ) ) 
C 

50  CONTINUE 
STOP 
END 
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IMSL  ROUTINE  NAME  >  ZSCNT 

PURPOSE  -  SOLVE  A  SYSTEM  OF  NONLINEAR  EQUATIONS 

USAGE  -  CALL  ZSCNT  (rCN,NSIG,N,ITMAX,PAR,X,FNORM, 

WK.IER) 

ARGUMENTS  FCN  -  THE  NAME  OF  A  USER-SUPPLIED  SUBROUTINE  WHICH 

EVALUATES  THE  SYSTEM  OF  EQUATIONS  TO  BE 
SOLVED.  FCN  MUST  BE  DECLARED  EXTERNAL  IN 
THE  CALLING  PROGRAM  AND  MUST  HAVE  THE 
FOLLOWING  FORM, 

SUBROUTINE  FCN (X,F,N,PAR) 

DIMENSION  X(N) ,F(N) ,PAR(1) 

Fd)- 


NSIG 

N 

ITMAX 

PAR 

X 

FNORM 

WK 

lER 


F(N)- 

RETURN 

END 

GIVEN  X(l) . . .X(N) ,  FCN  MUST  EVALUATE  THE 
FUNCTIONS  F(1)...F(N)  WHICH  ARE  TO  BE  MADE 
ZERO.  X  SHOULD  NOT  BE  ALTERED  BY  FCN.  THE 
PARAMETERS  IN  VECTOR  PAR  (SEE  ARGUMENT 
PAR  BELOW)  MAY  ALSO  BE  USED  IN  THE 
CALCULATION  OF  F(l)...F(N). 

THE  NUMBER  OF  DIGITS  OF  ACCURACY  DESIRED 
IN  THE  COMPUTED  ROOT  (INPUT). 

THE  NUMBER  OF  EQUATIONS  TO  BE  SOLVED  AND 
THE  NUMBER  OF  UNKNOWNS  (INPUT). 

THE  MAXIMUM  ALLOWABLE  NUMBER  OF  ITERATIONS 
(INPUT) . 

PAR  CONTAINS  A  PAR.AMrTER  SET  WHICH  IS 

PASSED  TO  THE  USER  SUPPLIED  FUNCTION  FCN. 
PAR  MAY  BE  USED  TO  PASS  ANY  AUXILIARY 
PARAMETERS  NECESSARY  FOR  CO.MPUTATION  OF 
THE  FUNCTION  FCN.  (INPUT) 

A  VECTOR  OF  LENGTH  N.  (INPUT/OUTPUT)  ON  INPU”^ 
X  IS  THE  INITIAL  APPROXIMATION  TO  THE  ROOT 
ON  OUTPUT,  X  IS  THE  BEST  APPROXIMATION  TO 
THE  ROOT  FOUND  SV  ZSCNT. 

ON  OUTPUT,  FNORM  IS  EQUAL  TO 

F(i) **2+. . .F (N) AT  THE  POINT  X. 

WORK  VECTOR  OF  LENGTH  (N-1)  *  ( 3*N>8) 

ERROR  PARAMETER.  (OUTPUT) 

TERMINAL  ERROR 

lER  «  129  INDICATES  THAT  ZSCNT  FAILED  TO 
CONVERGE  WITHIN  ITMAX  ITERATIONS.  THE 
USER  MAY  INCREASE  ITMAX  OR  TRY  A  NEW 
INITIAL  GUESS. 

lER  «  130  INDICATES  THE  ALGCRITHM  WAS 
UNABLE  TO  IMPROVE  ON  THE  RETURNED  VALUE 
OF  X.  THIS  SITUATION  ARISES  WHEN  THE 
SOLUTION  CANNOT  BE  DETERMINED  TO  NSIG 
DIGITS  DUE  TO  ERRORS  IN  THE  FUNCTION 
VALUES.  IT  MAY  ALSO  INDICATE  THAT  THE 
ROUTINE  IS  TRAPPED  IN  THE  AREA  OF  A 
LOCAL  MINIMUM.  THE  USER  MAY  TRY  A  NEW 
INITIAL  GUESS. 
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PRECIS lON/HAROWARE  -  SINGLE  AND  DOUBLE/H32 

-  SINGL£/H36«H48,H60 

REQD.  IMSL  ROUTINES  -  SINGLE/GGUBFS ,LEQT2F, LUDATN, LUELMN, LUREFN, 

UERSET , UERTST , UGETIO , ZSCNU 

-  DOUBLE/GGUBFS , LEQT2F , LUDATN , LUELMN , LUREFN , 

UERSET, UERTST , UGETIO , VXADD , VXMUL , VXSTO , 
ZSCNU 


NOTATION 


Aloorithm 


-  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


ZSCNT  is  based  on  an  algorithm  described  in  the  reference  sited  below 
At  each  step  of  the  iterative  process,  there  are  n+1  trial  solutions 

Multipliers,  a ^ , j»l , 2 , . . . , n+1  are  determined  by  solving 
the  linear  system: 


n+1 

E  •  0 

j»i 


•*•^1,2,  •  • » ,n 


:hen  the  new  trial  solution  is  defined  by 


*  •  E  • 

i»l 


See  reference: 


Wolfe,  Phillip,  "The  Secant  Method  for  Simultaneous  Non-Ii 
tions* ,  Communications  of  the  a£m,  Vol.  ~1',  rT5T5TT“ 


-linear  ioua- 


Example 


A  3  by  3  system  of  non-linear  equations  is  to  be  solved.  The  values 
of  the  parameters  in  PAR  are  passed  to  FCN  by  the  calling  program. 

Input: 

INTEGER  N,ITMAX,IER 

REAL  WK(68) ,X(3) ,PAR(3) ,FNORM 

DATA  PAR/27.,10.,7./ 

EXTERNAL  FCN 


ZSCNT- 2 
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N  -  3 
NSIG  -  5 
UMAX  »  100 
X(l)  -  4. 

X(2)  •  4. 

X(3)  -  4. 

CALL  2SCNT  (PCN, NSIG, N,  ITMAX, FAR, X,FNORM,vnC,  lER 


STOP 

END 

SUBROUTINE  FCN (X,F,N,PAR} 

REAL  X(3} ,F(3) ,PAR(3) 

F(1)-X{1)+EXP(X(1)-1.)**-(X(2)+X(3)  )**2-PAR(l) 
F(2)-EXP(X(2)-2.)/X(l)+X(3)**2-PAR(2) 
F(3)«X(3)+SIN(X(2)-2.)+X(2)  **2-PAR{3) 

RETURN 

END 

Output : 

FNORM  -  .00000 

X(l)  -  1.00001 
X(2)  •  2.00000 

X(3)  -  3.00000 
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IMSL  ROUTINE 

PURPOSE 

USAGE 

ARGUMENTS 


NAME  •  ZSPON 

-  SOLVE  A  SYSTEM  OF  NONLINEAR  EQUATIONS 

-  CALL  ZSPOW  (FCN,NSIG,N,ITMAX,PAR,X,rNORM, 

WX,I£R) 

FCN  -  THE  NAME  OF  A  USER-SUPPLIED  SUBROUTINE  WHICH 
EVALUATES  THE  SYSTEM  OF  EQUATIONS  TO  BE 
SOLVED.  FCN  MUST  BE  DECLARED  EXTERNAL  IN 
THE  CALLING  PROGRAM  AND  MUST  HAVE  THE 
FOLLOWING  FORM, 

SUBROUTINE  FCN(X,F,N,PAR) 

REAL  X(N)  ,F(N)  ,PAR(1) 

F(l)- 

F*(N)« 

RETURN 

END 

GIVEN  X(1)...X(N),  FCN  MUST  EVALUATE  THE 
FUNCTIONS  r(l)...F{N)  WHICH  ARE  TO  BE  MADE 
ZERO.  X  SHOULD  NOT  BE  ALTERED  BY  FCN.  THE 
PARAMETERS  IN  VECTOR  PAR  (SEE  ARGUMENT 
PAR  BELOW)  MAY  ALSO  BE  USED  IN  THE 
CALCULATION  OF  F(1)...F(N). 

NSIG  -  THE  NUMBER  OP  DIGITS  OP  ACCURACY  DESIRED 
IN  THE  COMPUTED  ROOT.  (INPUT) 

N  -  THE  NUMBER  OF  EQUATIONS  TO  BE  SOLVED  AND 

THE  NUMBER  OF  UNKNOWNS.  (INPUT) 

ITMAX  -  THE  MAXIMUM  ALLOWABLE  NUMBER  OP  ITERATIONS. 

(INPUT)  THE  MAXIMUM  NUMBER  OF  CALLS  TO  FCN 
IS  ITMAX* (N+1).  SUGGESTED  VALUE  •  200. 

PAR  -  PAR  CONTAINS  A  PARAMETER  SET  WHICH  IS 

PASSED  TO  THE  USER-SUPPLIED  FUNCTION  FCN. 
PAR  MAY  BE  USED  TO  PASS  ANY  AUXILIARY 
PARAMETERS  NECESSARY  FOR  COMPUTATION  OF 
THE  FUNCTION  FCN.  (INPUT) 

X  -  A  VECTOR  OF  LENGTH  N.  ( INPUT /OUTPUT)  ON  INPUT, 

X  IS  THE  INITIAL  APPROXIMATION  TO  THE  ROOT. 
ON  OUTPUT,  X  IS  THE  BEST  APPROXIMATION  TO 
THE  ROOT  FOUND  BY  ZSPOW. 

FNORM  -  ON  OUTPUT,  FNORM  IS  EQUAL  TO 

F(l)**2+. . .F(N)**2  AT  THE  POINT  X. 

WK  -  WORK  VECTOR  OF  LENGTH  N*(3*N-15)/2 
lER  -  ERROR  PARAMETER.  (OUTPUT*) 

TERMINAL  ERROR 

lER  ■  129  INDICATES  THAT  THE  NUMBER  OF 
CALLS  TO  FCN  HAS  EXCEEDED  IT.MAX*(N-:). 

THE  USER  MAY  TRY  A  NEW  INITIAL  GUESS. 
lER  -130  INDICATES  THAT  NSIG  IS  TOO 
LARGE.  NO  FURTHER  IMPROVE.MENT  IN  THE 
APPROXIMATE  SOLUTION  X  IS  POSSIBLE. 

THE  USER  SHOULD  DECREASE  NSIG. 
lER  «  131  INDICATES  THAT  THE  ITERATION 
HAS  NOT  MADE  GOOD  PROGRESS.  THE  USER 
MAY  TRY  A  NEW  INITIAL  GUESS. 
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PRECISION/HARDWARE 


SINGLE  AND  DOUBLE/R32 
SXNGLE/H36,H48,H60 


REQD.  IMSL  ROUTINES  -  SlNCLE/UERTST,UGETIO,  VBLA-SNR.M2  ,  ZSPWA, 

ZSPWB , ZSPWC , ZSPWD , ZSPWE , ZSPWF , ZSPWG 
-  DOUBLE/UERTST , UGETIO , VBLA-DNRM2 , ZSPWA , 

ZSPWB , ZSPWC , ZSPWD , ZSPWE , ZSPWF , ZSPWG 

NOTATION  -  INFORMATION  ON  SPECIAL  NOTATION  AND 

CONVENTIONS  IS  AVAILABLE  IN  THE  MANUAL 
INTRODUCTION  OR  THROUGH  IMSL  ROUTINE  UHELP 


Algorithm 

ZSPOW  i*  based  on  the  MINPACK  subroutine  HYBRDl,  which  uses  a 
modification  of  M.J.D.  Powell's  hybrid  algorithm.  This 
algorithm  is  a  variation  of  Newton's  meth^  which  uses  a 
finite  difference  approximation  to  the  Jacobian  and  takes 
precautions  to  avoid  large  step  sizes  or  increasing  residuals. 

See  reference: 


1.  Mor4,  Jorge;  Garbow,  Burton;  and  Hillstrom»  Kenneth,  User 
Guide  for  MINPACK- 1,  Argonne  National  Laboratory  Report 
ANL-80-f4,  Argonne,  Illinois,  August,  1980. 

Example 

A  3  by  3  system  of  non-linear  equations  is  to  be  solved. 

Input : 

INTEGER  N,NSIG,ITMAX,IER 
REAL  PAR{3) ,X(3) ,rNORM,WK(36) 

DATA  PAR/27.,10.,7./ 

EXTERNAL  FCN 
N  «  3 

NSIG  -  4 

ITMAX  •  100 

X(i)  -  4.0 

X(2}  «  4.0 

X{3)  »  4.0 

CALL  ZSPOW  (FCN ,  NSIG ,  N ,  ITMAX ,  PAR,X ,  FNORM ,  WK ,  lER) 


END 

SUBROUTINE  FCN  (X,F,N,PAR) 

INTEGER  N 

REAL  X(N)  ,F(N)  ,PAR(1) 

F(1)«X(1)-EXP(X(1)-1.0)*(X(2)*X(3)  )•  (X(2)^X{3)  )-PAR(l) 

F(2)-EXP{X(2)-2.0)/X(l)+X(3)*X(3)-PAR(2) 

r{3)«X(3)+SIN{X(2)-2.0)+X(2)*X(2)-?AR(3) 

RETURN 

END 
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